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The mechanism of transition in the wake of a thin flat 
plate placed parallel to a uniform flow 


By HIROSHI SATO anp KYOICHI KURIKI 


Aeronautical Research Institute, University of Tokyo, Japan 
(Received 10 February 1961) 


A study was made of the laminar-turbulent transition of a wake behind a thin 
flat plate which was placed parallel to a uniform flow at subsonic speeds. 
Experimental results on the nature of the velocity fluctuations have made it 
possible to classify the transition region into three subregions: the linear region, 
the non-linear region and the three-dimensional region. 

In the linear region there is found a sinusoidal velocity fluctuation which is 
antisymmetrical with respect to the centre-line of the wake. The frequency of 
fluctuation is proportional to the 3 power of the free-stream velocity, and the 
amplitude increases exponentially in the direction of flow. The behaviour of 
small disturbances in the linear region was investigated in detail by inducing 
velocity fluctuation with an external excitation—actually sound from a loud- 
speaker. Solutions of the equation of a small disturbance superposed on the 
laminar flow were obtained numerically and compared with the experimental 
results. The agreement between the two was satisfactory. 

When the amplitude of fluctuation exceeds a certain value, the growth rate 
deviates from being exponential due to non-linear effects. Although velocity 
fluctuations in the non-linear region are still sinusoidal and two-dimensional, the 
experimental results on the distributions of amplitude and phase indicate that 
the flow pattern may be described by the model of a double row of vortices. 
This configuration lasts until three-dimensional distortion takes place in the 
final subregion, the three-dimensional region, in which the fluctuation loses 
regularity and gradually develops into turbulence without being accompanied 
by abrupt breakdown or turbulent bursts. 


1. Introduction 


The laminar-turbulent transition of boundary layers is a well-known pheno- 
menon which has been investigated theoretically and experimentally since the 
beginning of this century. The detailed process of transition, however, still 
remains unclarified. In the initial period of investigation, interest lay mainly in 
determining the transition Reynolds number, which was important from the 
practical viewpoint. Modern developments in experimental techniques of flow 
measurement have made it possible to make detailed observations of velocity 
fluctuations in the boundary layer. Thus, investigations in the last decade have 
been focused on the understanding of the mechanism of transition. Recent 
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experimental results have indicated that two typical patterns of velocity 
fluctuation are essential in the transition. One is the regular sinusoidal fluctuation 
of small amplitude, and the other is the instantaneous turbulent burst of 
much higher intensity. 

Sinusoidal fluctuations were first discovered experimentally by Schubauer & 
Skramstad (1948) in the boundary layer along a flat plate. They obtained good 
agreement between their results and the linearized stability theory which had been 
developed by Tollmien, Schlichting and others. Since then, extended studies 
have been made on the properties of sinusoidal fluctuations both experimentally 
and theoretically. Sinusoidal velocity fluctuations have been found in various 
flow fields at subsonic and supersonic speeds, and agreement between linearized 
theory and the experimental results is satisfactory (Sato 1959a, b; Laufer & 
Vrebalovich 1960). It is now well established that the behaviour of fluctuations 
of small amplitude in the boundary layer is described by the linearized theory. 
On the other hand, our knowledge on the turbulent burst is poor. Experimental 
results obtained by various investigators sometimes show disagreements, and 
no concrete conclusions have been obtained on the origin and the law of develop- 
ment of a burst. Moreover, in the transition region of free boundary layers, no 
bursts have been observed so far (Sato 1959a, b). In order to obtain a full under- 
standing of the detailed mechanism of transition, more extensive studies are 
wanted. 

The present investigation has been undertaken with the intention of clarifying 
the process of transition in wakes. A thin flat plate was placed in a uniform flow 
with zero angle of attack so as to produce a simple, well-defined wake. It is 
obvious that the wake of a cylinder is inadequate for detailed studies of transi- 
tion since the flow field behind a cylinder is much too complicated owing 
to boundary-layer separation and rolling up of vortices in the wake. The 
velocity profile of a laminar wake was calculated by Goldstein (1933), and 
stability calculations based on the Orr-Sommerfeld equation were carried out by 
Hollingdale (1940) and McKoen (1956) in the case of large Reynolds number. 
It was Hollingdale (1940) who made the first systematic observation of velocity 
fluctuations in the wake of a flat plate by taking pictures in a water tank. 
Taneda (1958) has extended these observations to a wider range of Reynolds 
number with an improved experimental arrangement. He has found sinusoidal 
velocity fluctuations whose frequency is proportional to 3 power of the flow 
speed. The present work was started in order to make a more detailed survey of 
velocity fluctuations in wakes, using modern techniques of hot-wire anemo- 
metry in a low-turbulence wind-tunnel. The first part of this paper gives a 
general view of the transition phenomenon. The second part deals with the 
sinusoidal velocity fluctuations, and the third and fourth parts are devoted to 
clarifying the role of non-linear effects and three-dimensional distortion in the 
process of development into turbulence. 
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2. Experimental arrangement 


The experiment was conducted in a 60 x 60 cm low-turbulence wind-tunnel 
which was newly constructed at the Aeronautical Research Institute. The 
tunnel is of non-return type with the test-section downstream of the fan. The 
contraction ratio, i.e. the ratio of cross-sectional area of the settling chamber to 
that of the test-section, is 16:1. Eleven damping screens were installed in the 
settling chamber. These arrangements, incorporated with a proper design of air- 
intake, fan and diffuser, diminished the turbulence level at the test-section 
to 0-05°% at a flow velocity of 10 m/sec. Background noise and mechanical 
vibrations were carefully reduced by suitable construction, since they have a 
great influence on the transition process. The wind speed at the test-section was 


3m 
Y Wake | 
Plate 
30cm" Hot-wire 
anemometer 
Traversing mechanism 
Loudspeaker 


Ficure 1. Layout of test-section. 


variable from 1 to 30 m/sec. The general layout of the test-section is shown 
in figure 1, in which X is the distance downstream from the trailing edge of the 
plate, Y is measured perpendicular to the plate, and Z is the spanwise distance 
from the central point of the trailing edge. 

The side walls of the test-section were made of aluminum plates which were 
adjusted in order to keep the streamwise static-pressure distribution constant. 
One of the side walls had slits parallel and normal to the flow which allowed the 
insertion and traversing of hot-wire anemometers and pressure probes in the 
X-, Y- and Z-directions. 

The two-dimensional wake was realized behind a thin flat plate which was 
spanned vertically along the centre-line of the wind-tunnel. The parallel align- 
ment between the plate and the direction of flow was accomplished by measuring 
the mean-velocity distribution at the trailing edge of the plate. Three kinds of 
flat plate were used in order to clarify the effect of the thickness of the plate, 
as well as the sharpness of trailing edge, on the transition process. Dimensions 
of these plates are summarized in table 1. 

The chord and span of all plates were 30 and 60 cm, respectively. Plate I was 
made of aluminium, carefully machined into a thin airfoil over the whole length 
and polished to a mirror-like surface. Plate II was a thin brass foil without 
machining. Plate III] was made of aluminium and was sharpened at the leading 
edge, with the remaining part kept at the same thickness until the square 
trailing edge. The Reynolds number based on the length of plate ranged from 
6 x 104 to 4x 10°. 
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For the measurement of the mean-velocity distribution, both fine Pitot tubes 
and hot-wire anemometers were used, although the measurements by Pitot 
tubes formed a very small part of the present experiment. For the hot-wire 
anemometers, use was made of 10° rhodium-platinum wire, 5 ~ in diameter 
and about 2mm in length. For the measurement of the w-component of 
velocity, the hot-wire was rotated around the axis of support. 


Thickness 
Maximum Thickness at the 
thickness ratio trailing edge 
(mm) (%) (mm) 
Plate I 3 it 0-1 
Plate IT 0:3 0-1 0:3 
Plate III 3 1 3 


TABLE 1. Dimensions of plate. 


The mean velocity and the root-mean-square of the fluctuating velocity were 
measured by conventional hot-wire equipment which has been previously 
reported (Sato, Kobashi, Iuchi, Yamamoto & Onda 1954). The spectral distri- 
bution of the velocity fluctuation was observed by a band-pass filter. The mean- 
cube of the fluctuation was measured by a newly constructed unit in which the 
cubic curve was approximated by the superposition of the plate currents of ten 
triodes which were properly biased and loaded. A mean-square output meter 
was constructed on the same principle. This type of circuit had the advantages 
of higher input impedance and lower power consumption compared with similar 
circuits constructed by diodes. The fluctuation patterns were observed and 
recorded by a dual-beam cathode-ray oscilloscope. 

The anemometers were traversed in the X-, Y- and Z-directions. The X- and 
Z-positions were determined by scales, and the Y-position was indicated by a 
precision dial gauge. The accuracy of positioning was 0-2 mm in the X- and 
Z-directions, and 0-01 mm in the Y-direction. 

The flow was artificially excited by the sound radiated from a 10 W loud- 
speaker placed at the exit of the test-section as shown in figure 1. The maximum 
available intensity of the sound was about 100 db, although a much weaker 
sound was enough to induce transition in the laminar wake. The excitation 
frequency was varied from 500 to 1000 c/s, excluding the resonance frequency of 
the wind-tunnel, which was about 200 c/s. The intensity of sound was uniform 
throughout the test-section, and no resonance effect was detected in the above- 
mentioned frequency range. 


3. General aspects of transition 


At first, measurements were made in the boundary layer along the surface 
of the plate. The streamwise static-pressure variation along the surface was very 
small because the plates were very thin. The boundary layer was laminar every- 
where on the plate, and no traces of a velocity fluctuation were found even at the 
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highest wind-speed. The velocity profile was of Blasius-type, and the streamwise 
development of the boundary layer was in good agreement with the well known 
theoretical calculation for a flat plate. These results were quite satisfactory as 
regards fixing the initial condition of the laminar two-dimensional wake. 

After separation at the trailing edge, the boundary layers from both sides 
coalesced and formed a wake, which still remained laminar to some distance 
downstream. The velocity profile gradually changed from Blasius-type into 
laminar-wake type. Velocity distributions at various X-stations are shown in 
figure 2, where U, denotes the free-stream velocity. At X = 0, the distribution 
isof Blasius-type. Until X = 30 or 40 mm, the distribution varies slowly, while 
a sharp increase of central velocity is found from X = 40 to 60mm. In the 
figure two unusual facts are demonstrated. First the velocity at the outer 
part of the wake sometimes exceeds the free-stream velocity (for example, at 
X = 120mm, Y = 4to5 mm). Secondly, the central velocity is approximately 
equal at X = 60, 120 and 150 mm. 


A Ach. 


° X=30mm 
— —— 
X=60mm 
4 X=120mm 
4 X=150mm 


Y (mm) 


Figure 2. Mean-velocity distribution. Plate I, Ug = 10-0 m/sec. 


In order to illustrate the streamwise variation of velocity distribution more 
clearly, the velocity U,on the centre-line is plotted against X in figure 3. Curves 
for various experimental conditions show the same trend. The theoretical curve 
for the laminar wake given by Goldstein (1933) is in agreement with the experi- 
mental results when X is small. As X is increased, the experimental curves 
deviate from Goldstein’s result, reach maxima, decrease and again increase 
gradually. The point of deviation from the theoretical curve may be called the 
‘transition point’. It must be kept in mind, however, that the transition never 
takes place at a point and, moreover, that the definition of transition itself is 
arbitrary. A more detailed account of the velocity distribution in the transition 
region will be given in the following paragraphs. 

At X = 5 to 30 mm, depending on experimental conditions, there appears a 
velocity fluctuation whose wave-form is regular and sinusoidal. By careful 
investigation it was confirmed that the observed sinusoidal fluctuation was not 
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the result of either mechanical vibrations of the hot-wire or of acoustical 
resonance in the wind-tunnel. The frequency of fluctuation was the same 
throughout the whole region in which the wave-form was sinusoidal. The 
amplitude increased downstream until the wave-form was distorted into an 
irregular pattern. Since similar sinusoidal fluctuations have been already found 
in separated layers (Sato 1956, 19596) and jets (Wehrmann & Wille 1958; Sato 
1959a), the existence of sinusoidal fluctuations in the wake is not surprising. 
It has also been pointed out that the frequency of such fluctuations is propor- 
tional to the 3 power of the flow velocity according to dimensional reasoning 
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FicurRE 3. Streamwise variation of mean velocity U, on the centre-line 
under various experimental conditions. 


(Sato 1956). The logarithmic plot of frequency in the wake verified the same 
power law, as is shown in figure 4. The difference in frequency for three different 
plates was due to the difference in the thickness of the trailing edge. For a 
thick trailing edge, such as on Plate III, the thickness of the wake is com- 
paratively large, and therefore the frequency of fluctuation is lower. 

The development of the velocity fluctuations is illustrated in figures 5 and 6, 
which show results obtained with slight artificial excitation at the ‘natural’ 
frequency (i.e. the frequency of the sinusoidal fluctuations observed in the 
natural transition). Since the excitation was very slight, it had little effect on 
the transition process. On the other hand, the excitation improved the repro- 
ducibility of the flow field, and it was also used as a phase standard in the 
measurement of the phase relation of the velocity fluctuations. A detailed 
discussion concerning the interrelation between natural and excited transition 
is given in § 7. The wave-form of the u-fluctuation (figure 5) shows a gradual 
development from sinusoidal into irregular patterns. No turbulent bursts or 
momentary breakdowns have been observed. This type of transition is in 
contrast to that of a boundary layer on a solid wall, and is in accordance with 
transitions in separated layers and jets. Figure 6 is a map of fluctuation patterns 
at various X- and Y-positions. The sinusoidal fluctuation was observed from 
X = 20 to about 400mm. Near the centre-line the second-harmonic component 
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prevailed as shown in figure 5. The irregular fluctuation first appeared near the 
centre-line and spread, as X was increased, until about X = 400 mm, where 
the whole wake was covered by irregular fluctuations. 

It is worth while to note that at low Reynolds numbers the velocity fluctuation 
decays by viscous dissipation before the wave-form becomes irregular. In this 
case, no turbulent region exists in the wake even though the intensity of sinu- 
soidal fluctuation may be considerable. This is another sharp contrast to the 
transition of a boundary layer on a solid wall, in which the continuous energy 
supply due to Reynolds stresses and velocity gradient is maintained. In the low 
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Ficure 4. Frequency of natural sinusoidal fluctuations plotted 
against free-stream velocity Up. 


Reynolds number wake, the velocity gradient as well as the Reynolds stresses 
becomes progressively smaller downstream, and the viscous dissipation becomes 
pronounced. At higher Reynolds number, the wave-form of the fluctuations 
becomes irregular before it dies out, but it is still doubtful whether the wake 
becomes ‘fully turbulent’ or not. Moreover, even when the central part of the 
wake becomes turbulent, the regular fluctuation can persist in the outer part, as 
shown in figure 6. In these cases the definition of transition is vague. If the 
existence of a fully-turbulent region is a prerequisite of transition, there occurs 
no transition in the above-mentioned cases. Nevertheless, a high-level regular 
fluctuation exists in the wake and the mean-velocity distribution is different 
from that of laminar wake. We might call the region which is neither laminar nor 
turbulent the ‘transition region’. In so far as we use the term ‘transition’ in this 
sense, the transition point might be the end point of the laminar region and not 
the beginning point of turbulent region. 

The intensity of the w-fluctuations, with Plate I and U, = 10-0 m/sec, is 
plotted in figure 7. When X was small, the distribution was almost similar, and 
positions of peak values approximately coincided with positions of maximum 
gradient in the mean-velocity distribution. As X was increased, the intensity 
increased until about X = 60mm and then decreased. Distributions from 
X = 60 to 150 mm have maxima on the centre-line. This is due to the increase 
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of the second-harmonic component, which is shown in figures 5 and 6. The 
intensity of fluctuation in the present experiment was rather small compared 
with that in the wake behind a cylinder. This fact might be explained by the 
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Ficure 5. For legend see facing page. 
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X=300mm, Y=1mm 


X=300 mm, Y=5mm 


X=500mm, Y=3mm 
FicurE 5. Oscillographie records of w-fluctuations in wake. Plate I, Ug = 10-0 m/sec. 
Time goes from left to right and time interval between dots is 0-01 sec. 


- 8mm 


=> Sinusoidal wave we 4 
Laminar 
0 
Trailing edge \ 300mm 
= 00 
-8 
Ficure 6. Map of fluctuation patterns. Plate I, Uy = 10-0 m/sec. 


disturbance at the trailing edge of a flat plate being smaller than that behind a 
cylinder. 

After these investigations were made on the velocity fluctuation in the 
transition region, we found it convenient to consider a division into three sub- 
regions, i.e. the linear region, the non-linear region, and the three-dimensional 
region, which are characterized by the pattern of fluctuations. 

1. Linear region. In this region were observed two-dimensional sinusoidal 
velocity fluctuations which were amplified exponentially downstream. Referring 
to the previous investigations on free boundary layers (Sato 1959a, b), it is 
expected that these sinusoidal fluctuations may be described by the linearized 
theory. This region ends at the X-station where the amplitude of fluctuation 
reaches the ‘linearity limit’: in other words, where the law of amplification 
deviates from being exponential. The mean-velocity distribution in this region 
is in agreement with Goldstein’s theoretical result. 

2. Non-linear region. As a result of amplification of the sinusoidal fluctuations, 
second and higher harmonics are found in the wave-form, as shown in figure 5. 
Two-dimensionality is still maintained, that is, w? is small compared with w? 
and v? and the phase difference of the velocity fluctuations in the Z-direction is 
small. The variation of mean velocity is remarkable in this region, being the 
result of intense energy supply from the mean flow to the velocity fluctuations. 

3. Three-dimensional region. Further downstream the velocity fluctuation 
becomes three-dimensional, while some periodicity still persists. The wave-form 
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is irregular in the central part and is sinusoidal at the outer part of the wake. 
The turbulent region starts where the remaining regularity diminishes. 

In the following paragraphs the mechanism of transition will be made clearer 
by describing in detail these three regions. 
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Figure 7. Distribution of intensity of w-fluctuations. Plate I, Uy = 10-0 m/sec. 


4. The linear region 


As shown in the previous paragraph, sinusoidal fluctuations are found in the 
transition region of the wake. They are weak at first, grow as the distance down- 
stream is increased, and finally develop into irregular fluctuations still further 
downstream. The nature of the sinusoidal fluctuations has been clarified by both 
experimental and theoretical investigations. 


4.1. Experimental results 


Experiments on sinusoidal fluctuations have been made in two cases, one for 
natural transition and the other for transition caused by the external artificial 
disturbance. The latter is important for the experimental verification of 
linearized theory, which in turn is useful for clarifying the mechanism of the 
linear region in natural transition. Since natural transition results from the 
development of small unknown and uncontrollable disturbances, there is no 
essential difference between the natural and excited transition except the nature 
of the original disturbances. However, before detailed measurements are made 
with artificial excitations, it should be ascertained that the flow field is not 
changed too much by them. If the intensity of the artificial disturbance is too 
high, non-linear effects set in at X-stations that would lie in the linear region of 
natural transition. On the other hand, if the incoming sound is too weak, the 
induced fluctuation is masked by natural fluctuations, and in this case accurate 
measurements are difficult. The proper level of exciting sound was carefully 
decided as a compromise between these two considerations. 
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Experiments were carried out using Plate I at U, = 10-0 m/sec. The frequency 
of natural fluctuations in this case was around 730 c/s. The frequency range of 
the artificial excitation was from 480 to 850 c/s. Outside this range the observed 
fluctuation was hardly affected by the sound. 

The phase of the sinusoidal fluctuations was not constant in the X-direction, 
that is, the fluctuation was in fact a travelling wave. The propagation velocity 
was measured by moving a hot-wire anemometer in the X-direction. Since the 
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FicuRE 8. Propagation velocity of sinusoidal waves induced by artificial excitation. 
Plate I, Up = 10-0 m/sec. 


variation of phase in the Y-direction was remarkably large near the centre-line, 
as shown later, the phase difference in the X-direction had to be determined at 
large Y. Results for various frequencies of excitation are shown in figure 8, in 
which the propagation velocity c,, expressed as a fraction of U,, is plotted 
against X. The value of c, increases as X is increased, reaches a maximum at 
about X = 60 to 70 mm, and settles at about 0-9 for larger X. It is interesting to 
note that the variation of c, in the X-direction is similar to that of U,. Detailed 
discussion of the values of ¢, is found in §4.3. 

The increase in amplitude of the spectral components in the X-direction is 
illustrated in figure 9, in which (u?)b ax is the maximum value of the Y-distribu- 
tion, of root-mean-square intensity of the spectral component at each X-station 
in the presence of excitation, whose frequency was as indicated. The value was 
obtained by adjusting the central frequency of a band-pass filter to be exactly 
equal to the frequency of the sound excitation, and then traversing the hot-wire 
in the Y-direction. Since (u?)b ax is plotted on an arbitrary scale, comparison 
between values at different frequencies would be meaningless, although the 
relative values for the same frequency are correctly expressed. Each spectral 
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component was amplified exponentially at small X and deviated from the 
exponential form at X = 30 to 35mm due to non-linear effects. From the 
gradient of each curve the rate of amplification in the X-direction is calculated. 
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FicurE 9. Growth of spectral components. Plate I, Uy) = 10-0 m/sec. 


Distributions in the Y-direction of amplitude and phase of natural sinusoidal 
fluctuations are shown in figure 10. The amplitude distribution has two peaks at 
about Y = +1 mm, and becomes nearly zero on the centre-line. A change in 
phase is observed at -3 mm < Y < 3mm. This is the reason why the variation 
of phase in the X-direction had to be measured at Y = 4or5 mm. Although a 
phase shift of about 70 degrees existed between Y = —1 and —5 mm, the 
phase difference at any two symmetrical points with respect to the centre- 
line was always approximately 180 degrees. These results are to be compared 
presently with the linearized theory. 

Figure 11 shows distributions of amplitude and phase of the second-harmonic 
component (1460 c/s) which is found predominantly near the centre-line, as 
described before. The amplitude is a maximum at Y = 0, and the phase distri- 
bution shows that this component is symmetrical with respect to the centre- 
line. For the generation of the harmonic component the non-linear effect might 
be responsible. A full discussion of this point is given in the following sections. 


4.2. Theoretical considerations 
A theoretical investigation of the stability of wakes to small disturbances has 
been carried out by Hollingdale (1940) and McKoen (1956). Their work is 
restricted to the case of neutrally stable disturbances and, moreover, their 
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the methods of approximation are rather crude. Therefore, in order to compare the 
the foregoing experimental results with theory, calculations extending to the non- 

ted. neutral case have been made with improved accuracy. Only the case of infinite 


Reynolds number has been considered. 
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FicurE 10. Amplitude and phase distributions of fundamental component. Plate I, 
U, = 10-0 m/see, X = 40 mm. 
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We follow the procedure of conventional linearized stability calculations. The 
1as basic flow is assumed to be two-dimensional and a function of y only, and lateral 
is disturbance velocities are neglected. The disturbance velocities u and v are 
eir assumed to be functions of x, y and time ¢. Quantities in the equations of 
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motion and continuity are made non-dimensional by taking as units of velocity 
and length the free-stream velocity U) and the half-breadth 6 of the mean- 
velocity distribution. If X, Y, Z and T are the dimensional coordinates and 
time, the corresponding dimensionless quantities are written x = X/b,y = Y/b, 
z= Z/bandt = TU,/b. 

We define astream-function y such that the velocity perturbations are given by 


and we assume that w and v are small enough to justify linearization of the 
equations of motion. We further assume that the Reynolds number is large 
enough for the viscous terms in the equations of motion to be neglected. (The 
justification for this is based on the fact that the stability characteristics of free 
boundary layers become insensitive to viscous effects when the Reynolds number 
is fairly high.) Then, assuming y to take the form 


= Py) exp (1) 


in which a is the wave-number and ¢ = c,+ic;, where c, is the propagation 
velocity and c; a measure of the rate of amplification of the disturbance, we find 
that ¢ satisfies the equation 


U 


which is the ‘inviscid’ form of the well-known Orr-Sommerfeld equation. 
The boundary conditions are 


2(0)=1, (0) = (3) 
+00) = 0, 


and the mean velocity distribution is assumed to be 


— ary 


in which a is taken as 0-69315 so as to make (U,—U)/(Uj)—U,) = } at y = +1. 
The validity of this expression is indicated in figure 13. Since, in a wake, sym- 
metrical disturbances are more stable than antisymmetrical ones, calculations 
have been made only for the latter. 

In order to find a set of eigenvalues which satisfies the boundary conditions, 
numerical integrations have been repeated by a high-speed electronic computer. 
The computed results for the eigenvalues a and c, taking (U,—U,)/U, = 0-692, 
are shown in figure 12. 

Before we compare these theoretical results with experiment, we have to 
make a transformation of time axis and X-axis according to the relation a = ¢,t, 
since the flow in a wind-tunnel is stationary and the disturbance is amplified as 
it travels in the X-direction. Therefore the spatial rate of amplification ac;/c, 
must be used instead of ac; for the comparison with experimental results. 
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4.3. Comparison of experimental results with linearized theory 


A non-dimensional plot of the velocity defect in the linear region of the wake is 
shown in figure 13. The theoretical distribution for a fully developed laminar 
wake is in the form exp (—ay?), which is indicated by a solid line. The velocity 
distribution at X = 30 mm according to the calculation by Goldstein (1933) is 
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Ficure 13. Non-dimensional mean-velocity distribution in linear region. 


as shown by a broken line. The two theoretical curves are very close, and the 
experimental data are in good agreement with these curves. Since the distribu- 
tion given by the broken line has no simple analytical expression, the exponential 
distribution is preferable as a basic flow for the calculation of disturbance 
properties. 

The non-dimensional frequency ac, of the observed natural fluctuations is 
plotted against Reynolds number in figure 14. Since the theoretical value of 
ac, for neutral oscillations at infinite Reynolds number is 0-92, the frequency of 
the natural fluctuations can be assumed to lie inside the unstable zone of the 
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theoretical stability diagram, although the curve of neutral stability for finite 
Reynolds numbers has not been derived. To the right of the figure is shown the 
theoretical spatial amplification rate at infinite Reynolds number, which has 
a maximum at about ac, = 0-5. Thus it is confirmed that the frequency of the 
sinusoidal fluctuations found in natural transition approximately corresponds 
to the frequency of the small disturbance which receives maximum spatial 
amplification according to the linearized theory. 


1:0 1-0 
ac, 
Neutral 
(R=) 
05 05 
0 0 J 
0 500 1000 0 0-2 
R= ac;/c, 


Figure 14. Comparison of theoretical and experimental results. Frequency of velocity 
fluctuation, ac, = 2mfb/Up. 


When sound waves were used as artificial excitation, it was impossible to 
observe the damping of disturbances just beyond the critical condition since the 
induced velocity fluctuations were not sufficiently strong. Therefore the neutral 
point was determined experimentally by the method introduced in a previous 
paper (Sato 19596). The ratio of the intensities of spectral components with and 
without excitation is plotted against non-dimensional frequency in figure 15. 
The neutral frequency is determined as the point where the ratio becomes unity. 
Arrows on the abscissa indicate the theoretical values which have been obtained 
by various investigators: e.g. by Hollingdale (1940) who obtained two values with 
different methods of approximation, and by McKoen (1956). The present value 
of ac, determined by numerical integration is the smallest among the four given 
values. The experimental results obtained under four different conditions are 
fairly consistent, and they agree most closely with the present theoretical 
estimate. 

Figure 16 shows the propagation velocity and spatial rate of amplification 
which are obtained from experimental data such as figures 8 and 9. The solid 
lines are theoretical results for amplified disturbances at infinite Reynolds 
number. The agreement between theory and experiment is fairly good in both 
c, and ac,/c,. 

Distributions of amplitude and phase in Y-direction are compared in figure 17. 
Experimental points denote the root-mean-square values of spectral components 
of the u-fluctuations in the natural transition, taking the maximum value as 
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ite unity and, for the phase angle of the spectral components, taking the value at 
the Y =—5 mm as zero. Theoretical results are given for amplified disturbances for 
has which 4(y) is complex. The amplitude function and phase angle for the u-fluctua- 
the tions were calculated from the real and imaginary parts of ¢’(y). As shown by a 
nds solid line in the figure, a gradual phase shift in the y-direction exists for the 
tial 
{ 
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sity Figure 15. Determination of frequency of neutral oscillation. Plate T. and 
denote root-mean-square values of spectral components of frequency f with and without 
artificial excitation respectively. 
to FicureE 16. Propagation velocity c, and spatial rate of amplification ac;/c,. 
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Ids amplified disturbance, in contrast to the case of neutral disturbances. Theoretical 
oth eigenvalues in the present case are a = 0-832, c; = 0-168, c, = 0-692, which 
correspond to a disturbance with approximately the maximum spatial rate of 
17. amplification. The agreement between theory and experiment is excellent. 
nts Thus the experimental results on the properties of sinusoidal fluctuations in 
aad the linear region are properly described by the linearized theory. Considering 
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the previous results on the separated layer and jet, it is now obvious that in the 
transition region of a free boundary layer there exists a linear region in which 
sinusoidal fluctuations grow as predicted by the theory of the Orr-Sommerfeld 
equation. 


5. The non-linear region 


The upstream boundary of the non-linear region is the X-station where the 
growth rate of velocity fluctuations deviates from being exponential. This 
deviation demonstrates the role of non-linear terms in the equation of motion. 
This initial point of the non-linear region might be called the ‘transition point’ 
since the non-linearity is an essential feature of turbulent motion. On the other 
hand, the transition point defined as the end-point of laminar region in the 
preceding section is determined from, for instance, the streamwise variation of 
velocity on the centre-line such as is shown in figure 3. These two kinds of transi- 
tion point, one defined from the velocity fluctuation and the other from the 
time-mean velocity, are found experimentally to be very close. This fact has 
been already pointed out in the case of separated layers and jets (Sato 1956, 
1959a). From the experimental results it has been confirmed that the two- 
dimensionality of the fluctuations is still maintained in the non-linear region. The 
downstream boundary of the non-linear region is defined as the X-station where 
the fluctuations become three-dimensional so that, for instance, u? and w? become 
comparable. With Plate I and U) = 10-0 m/sec, the extent of the non-linear 
region was found to be about X = 40 to 150 mm. 


5.1. Non-linear development of sinusoidal fluctuations 


As shown in figure 5, the non-linear effect on the wave-form of the fluctuations is 
observed in the vicinity of the centre-line as an increase of low-frequency irregular 
fluctuations. In the outer part of the wake the sinusoidal fluctuation persists with 
the same frequency. The non-linear effect seems to be slight between X = 40 
and 150mm as far as the wave-form is concerned. However, remarkable 
changes take place in the distributions of amplitude and phase of the spectral 
components, as is shown in figures 18 and 19. In figure 18, the root-mean-square 
values of the spectral component at the fundamental frequency are plotted on an 
arbitrary scale, while the relative values at various X-stations are correctly 
expressed. The distribution changes in a rather peculiar way as X is increased. 
From X = 40 to 80mm, the variation is gradual. The positions of the two 
peaks move outwards, maintaining a roughly similar type of distribution. At 
X = 120 and 150mm, the distribution shows three minima, at Y = 0 and 
Y =+25mm. Referring to figure 19, we see that a 180-degree phase reversal 
takes place at these three points. At each X-station, 180-degree reversal is 
observed at Y = 0, except at X = 80 mm which is in some sense a dividing line. 
Now it seems possible to classify the non-linear region into two parts. The 
upstream part more or less resembles the linear region, although the exponential 
amplification of fluctuations is no longer observed. In the downstream part, the 
distribution functions of amplitude and phase are completely different from 
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those in the upstream part. It is obvious that in the non-linear region a new 
equilibrium is established at about X = 100 mm. 

The next step of the investigation is to clarify the nature of the harmonic 
components in the non-linear region. An increase of harmonic components is 
expected as a result of non-linear interaction. In figure 20, the distribution of 
second harmonics is shown. While a remarkable increase in (w3,) takes place 
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Figure 18. Amplitude distribution of fundamental component in non-linear region. 
The ordinate scale is arbitrary. 
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Figure 19. Phase distribution of fundamental component in non-linear region. 


from X = 40 to 60mm, it decreases downstream from X = 60mm. This 
is a curious situation since the fundamental component still increases or 
only slightly decreases from X = 60 to 100mm. The ratio of harmonic com- 
ponent to fundamental component becomes a maximum at X = 60mm and 
decreases as X is increased further. Although non-linear effects might always be 
expected to increase the harmonic content, we here observe a contrary result. 
The distribution functic of the second-harmonic component in the 
Y-direction shows a reduction in the extent of the distribution from X = 60 
to 120 mm. This fact is also unusual. A’ more detailed discussion of figure 20 
is hardly justified, however, since the accuracy of measurement of the 
harmonic components was not high, owing to the unavoidable non-linearity of 
the characteristics of the hot-wire anemometer itself. The necessary correction 
for this non-linearity reached 50°, of the value of (w3,), while the harmonic 
component was very small compared with the fundamental component. 
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Another quantity which characterizes the non-linear distortion of sinusoidal 
fluctuations is (v3)'. When the wave-form remains symmetrical with respect to 
the zero line, the value of (w°)! is zero. Thus (wz), which was determir: 7 by the 
use of a cubing unit, is a measure of asymmetry of the wave-form. Due to the 
inevitable non-linearity of the hot-wire, as mentioned before, the measured value 
of (u?)} was not zero even for sinusoidal fluctuations when the amplitude was 
fairly large. Moreover, the wave-form distortion was usually accompanied by 
slow irregular velocity fluctuation. For these reasons, a high accuracy was not 
expected in the measurement of (w3)}. Figure 21 is the result obtained in the 
non-linear region. At X = 60mm, the distribution is different from those at 
X = 80 and 120 mm, which again suggests the peculiar character of the region 
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Figure 20. Amplitude distribution of second harmonies in the non-linear region. 
The ordinate scale is arbitrary. 


X = 80 to 150mm. [t is worth while to note that in the distribution at 
X = 120mn, the zeros of (w3)! at about Y = +2 mm coincide approximately 
with the phase-reversal points for the fundamental component. Outside 
Y =+2mm, the second-harmonic component is very small, so that non- 
linearity is manifested as the asymmetry of the wave-form. Inside Y = +2 mm, 
(w3,)} is nearly constant, whereas (u3)} changes sign. From these facts it is 
expected that the generation of harmonics and asymmetrical distortion of wave- 
form have no simple relation, although both are undoubtedly results of non- 
linear interactions. 

In order to clarify the situation at X = 120 mm, the non-dimensionalized 
mean-velocity distribution (figure 22) may be closely examined. We find two 
facts. One is the over-shoot of velocity between y = 1-6 and 2-5, and the 
other is the flat distribution near the centre-line. The distribution at X = 40 mm 
is added for comparison. The peculiarity of the distribution at X = 120 mm is 
clearly seen. 

5.2. Vortex model 


The curious features of the flow field in the region X = 40 to 150 mm are 
summarized as follows: 

1. Before the experiment was made, it seemed probable that in the wake the 
central velocity U, increases monotonically and approaches U, at large X. But 
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at about X = 80 mm, U, becomes a maximum, begins to decrease, and reaches 
a minimum at X = 120 mm before it gradually increases again. 

2. The mean-velocity distribution at X = 120mm is very different from 
that at smaller or larger X. The distribution has a flat part near the centre-line, 
and the velocity over-shoots the free-stream velocity at the edge of the wake. 

3. The fundamental component of the velocity fluctuations at X =-120 to 
150 mm becomes zero not only at Y = 0 but also at Y = +2 mm. The phase 
reversal also takes place there. 
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Ficure 21. Distribution of (%)) in non-linear region. 
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Figure 22. Non-dimensional mean-velocity distributions in linear and non-linear regions. 


4. As X is increased, the second-harmonic component increases until 
X = 60 mm, decreases, reaches a minimum, and increases a little again as X is 
increased further. The distribution in the Y-direction shows a systematic 
change over the range X = 40 to 150 mm. 

These features might be explained by an analytical method taking account of 
non-linear terms in the equation of motion. Since the linearized theory has had 
substantial experimental verification in the linear region, the mechanism of the 
non-linear region might also be clarified by an appropriate theoretical model. 
Practically, however, analysis of the non-linear equations of motion is very 
difficult, and no satisfactory theory has been reported so far. Thus we are 
inclined to take a tentative empirical model which explains the above-men- 
tioned experimental facts. 

In the past, many investigators have observed the formation of vortices in 
the wake of two-dimensional bodies, including a thin flat-plate (Hollingdale 
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1940; Taneda 1958). It is usually accepted that the sinusoidal fluctuations in 
the signal from a hot-wire anemometer corresponds to the moving vortex which 
is observed visually. In order to make such a correspondence apply in the 
present case, the physical properties posed on the vortex must be sufficiently 
complex to explain the wave-form of the hot-wire output. The usual simple 
vortex model is obviously inadequate. In order to approximate the linear 
region, a single row of vortices which travels in the X-direction as shown in 
figure 23(a) is useful. The adjacent vortices rotate in opposite direction and the 
radial distribution of circumferential velocity is assumed as shown. This model 
roughly explains the calculated and observed distributions of amplitude and 
phase of the velocity fluctuations in the linear region. The nature of the fluctua- 
tions in the non-linear region is explained as follows by the development of the 
single row into a double row as shown in figure 23(b). 


FicuRE 23. Vortex model. (a) Single row; (b) double row of vortices. 


In the single row, the superposition of vortices has no effect on the mean- 
velocity field. On the other hand, in the double-row configuration a velocity is 
induced in the negative X-direction between rows of vortices, as shown by the 
arrow. This results in a decrease of U,. Furthermore, a positive velocity is 
induced at the outer part of the wake. These facts explain the mean-velocity 
distribution at about X = 120 mm as shown in figure 22. In the double row of 
vortices, the phase inversion of the fundamental component takes place not only 
on the centre-line but on lines connecting centres of vertices in each row if the 
lateral spacing is large. The generation of second harmonics in the vicinity of 
the centre-line is also easily understood. Typical wave-forms of velocity fluctua- 
tions as expected from the vortex model are shown at the bottom of figure 23, 
depending on the location of hot-wire. These are in accord with the experimental 
observations. The change of sign of (wv)! on the line connecting centres of 
vortices is understood from these wave-forms. Thus the process in the non- 
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linear region is explained by the development of the single row into a double 
row of vortices, and if once the double row is established, the configuration is 
stable and remains unchanged until three-dimensional deformation begins. 

It is worth while to calculate an important geometrical parameter of a double 
row of vortices, that is, the ratio of lateral to longitudinal spacings. The value is 
0-281 in von Karman’s analysis, whereas in the present experiment it is-about 
0-39. This comparatively large lateral spacing allows the phase inversion of the 
fundamental component on the lines connecting centres of vortices. This fact is 
in contrast to the case of the wake of a circular cylinder, in which the phase 
reversal is found only at Y = 0 (Roshko 1954). 


5.3. Interaction between spectral components 
In the linear region, spectral components of different frequencies are damped or 
amplified independently. In the non-linear region, however, two spectral com- 
ponents may interact with each other through the non-linear terms in the 


© 730 c/s 
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0 05 1-0 
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Ficure 24. Variation of intensity of natural and excited spectral components in 
the non-linear region. X = 120mm. 


equations of motion. This interaction was experimentally investigated by 
superposing an artificial disturbance, that is, sound from the loudspeaker 
at a frequency other than that of the natural fluctuations. 

At first, as the result of non-linear interaction, we expected the generation of 
spectral components of frequency f, +f, where f, and f, denote the frequencies 
of the natural and artificial disturbances respectively. However, no peaks were 
observed experimentally at f,+/, in the energy spectrum of the u-fluctuations. 
By varying the intensity of sound from the loudspeaker, the wave-form and 
intensity of the velocity fluctuations in the non-linear region have been investi- 
gated. Figure 24 shows the change in the root-mean-square values of two 
spectral components. The induced spectral component (u?,)3 (590 c/s) increases 
monotonically as the input voltage to the loudspeaker is increased. On the 
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other hand, the component of natural fluctuation (730 c/s) diminishes when the 
intensity of excitation is increased. This ‘suppression effect’ is more remarkable 
when /;, is close to fy. This effect means, in other words, the supply of energy from 
the component at f, to that at f,. The detailed mechanism of the energy exchange 
is not yet known. 


(a) Without excitation 


(6) Input to loud-speaker 1 V. 


(c) Input to loud-speaker 2 V. 


(d) Input to loud-speaker 3 V. 
Figure 25. Oscillographic records of u-fluctuations in the non-linear region at X = 120 mm 
and Y = 5mm. Excitation 590 c/s. Upper trace is the output from the hot-wire anemo- 
meter, and lower trace is the input to the loudspeaker. Time goes from left to right and 
the time interval between dots is 0-01 sec. 


The wave-form in the presence of artificial excitation is shown in figure 25. 
Without artificial excitation the frequency of fluctuation is about 730 ¢/s (a). 
When the incoming sound is weak, the component of /, appears in a small 
fraction of time, while the f, component prevails (b). By increasing the intensity 
of sound, the probability of finding the f, component increases (c), and with the 
maximum intensity of sound the f, component is hardly found (d). From these 
records, it is obvious that the probability of finding induced fluctuation is a 
function of intensity of the artificial disturbance. No periodicity was found in 
the occasional change of the frequency of fluctuation. 


6. Three-dimensional region 

In the linear and non-linear regions, velocity fluctuations are approximately 
two-dimensional. The distortion of two-dimensional waves—in other words, the 
three-dimensional deformation of the row of vortices—takes place in the final 
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subregion of transition, the three-dimensional region, in which regularities in the 
fluctuations gradually diminish until the fully turbulent wake is formed. If the 
Reynolds number of flow is small, the fluctuation energy decays by viscous 
dissipation before the development into turbulence is completed. In this case, 
there is no turbulent region and the three-dimensional region is also hardly 
observed. 

With the use of Plate I at U, = 10-0 m/sec, the physical process which takes 
place in the three-dimensional region was as follows. 

The double row of vortices formed in the non-linear region became unstable 
at about X = 150mm. As a first sign of the deformation of the vortices, 
irregular wave-forms were observed near the centre-line of the wake while the 
fluctuation was still sinusoidal in the outer part. As X was increased, the region 
of irregular fluctuations extended into the outer part until the whole breadth of 
the wake was covered. At the same time, the phase relation of the sinusoidal 
fluctuations with respect to the Z-direction became irregular, and w-fluctuations 
became comparable in magnitude with u-fluctuation. The turbulent wake was 
thus established at about X = 400 mm, although the intensity of fluctuation was 
much less there compared with that in the non-linear region. No instantaneous 
turbulent bursts have been observed so far in this region. The change of wave- 
form was gradual. 

Before we zo into the quantitative discussion, we may suitably introduce some 
physical quantities to characterize the three-dimensionality. From the con- 
sideration that the difference in mean-velocity distribution at various Z-stations 
is small, three quantities concerning the velocity fluctuation have been used. 
They are the relative magnitude of (w?)} and (w#?)}, the phase relation in the 
Z-direction, and the distribution of spectral components at various Z-stations. 

The ratio (w?)}/(u2)} is plotted against Y in figure 26. At X = 150 mm, u- 
and w-fluctuations are nearly the same in magnitude, or (w?)! is even greater in 
the region —2mm < Y < 2mm, and outside this region (w)! is very small 
compared with (w?)!. Therefore, the irregular fluctuation observed in the central 
region is three-dimensional, but the sinusoidal fluctuation in the outer part of 
the wake is still two-dimensional. At X = 300 mm, the ratio is close to unity over 
the whole breadth of the wake. This fact resulted from the faster decay of (w?)! 
rather than the growth of (w?)! in the three-dimensional region. The value of 
(w?)4/(u?)? at X = 300 mm is close to that in the fully developed wake behind 
a cylinder. 

The straightness of the front of the sinusoidal wave is expressed by the phase 
relation in the Z-direction at large Y. The experimental result is shown in 
figure 27, in which the phase angle at Z = 0 is taken to be zero. At small X, the 
variation in the Z-direction is small. This means that the wave is approximately 
two-dimensional and propagates in the X-direction. At larger X, a phase 
difference of more than 180 degrees is observed at some Z-station, while the 
direction of propagation as a whole still remains in the X-direction. As mentioned 
before, the wave-form at large X is irregular near the centre-line of the wake. 
This means that the double row of vortices is first distorted in the region between 
rows. For instance, the phase of the velocity fluctuations at X = 250 mm could 


the 

ible 
rom 
nge 
3 
\ 

/\ 

mm 
mo- 
and 
25. 
(a), 

vall 
sity 
the 
ese 
Sa 
| in = 

ely 


346 Hiroshi Sato and Kyoichi Kuriki 


not be determined experimentally because of the irregularity in the central part 
of the wake. The phase distribution in the Z-direction at large Y might still be 
significant, since it expresses the three-dimensional distortion in a large scale. 
As shown in figure 27, the Z-station of the phase lag or advance is approximately 
fixed at various X-stations. This fact might have arisen from imperfections of 
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FicureE 26. Distribution of (w?)!/(@)? at X = 150 and 300 mm. 


90° 
X=60mm 
6 
9 X=150mm 


Z (mm) 


FicurE 27. Spanwise distribution of phase of fundamental component at various 
X-stations in the outer part of the wake. Arrow indicates the station where detailed 
measurements were made. 


the flat plate. Upon careful inspection of the plate, a slight waviness (about 
0-2 mm in height) was found in the trailing edge at about Z = —40 mm. Since 
all measurements in the linear and non-linear regions were made at Z = 10mm, 
as indicated by an arrow in figure 27, the experimental results described in 
preceding sections were uninfluenced by disturbances produced by this waviness. 
However, since the slightest spanwise irregularity can have a decisive effect on 
the breakdown of the laminar boundary layer on a solid wall (Schubauer 1958; 
Tani 1960), a measurement has been made on the distribution of spectral com- 


pol 
sta 
sta 
wit 
On 
4 | | col 
X=300mm | pal 
| d the 
far 
| 
1. 
mt 
7 ini 
th 
Ar 
| —100 —50 
| 100 
na 
an 
th 
| re 
7 re 
ca 
tu 
of 
di 
| H 
al 
re 


part 
ll be 
ale. 
tely 
s of 


Transition in the wake of a thin flat plate 347 


ponents at three Z-stations A, B and C as denoted in figure 27. A and C are 
stations of maximum lag and advance respectively, and B is the midway 
station. The experimental results obtained at A and C are in good accordance 
with those at Z = 10 mm where all systematic measurements have been made. 
On the other hand, the distribution at B with X = 150mm is somewhat 
different. The 180-degree phase shift is observed only at Y = 0 and the maxi- 
mum value of amplitude is about 1-3 times that at Z = 10mm. The phase 
configuration suggests the single-row configuration at B, although it is accom- 
panied by an appreciable distortion. At X = 200 mm, the difference between 
the three Z-stations is not observed any more. From the observations made so 
far, the peculiar region around station B has no decisive effect on the transition 
process in the whole wake. This fact is in contrast to the case of transition of the 
boundary layer along a solid wall, which depends to an important extent on the 
three-dimensionality in the layer. 

The three-dimensional region is succeeded by the turbulent region, on which no 
detailed measurements have been made in the course of the present investigation. 


7. Discussion 

First of all, the relation between natural and artificially excited transition 
must be clarified, because in the present experiment the transition was often 
initiated by sound from a loudspeaker. Generally speaking, in wind-tunnels 
there are various kinds of residual disturbance such as the free-stream turbu- 
lence, acoustic noise, mechanical vibration, ete., which undoubtedly have great 
influence on the transition process. These residual disturbances usually change 
day by day, and the reproducibility of transition pattern is sometimes very poor. 
And, moreover, the magnitude and nature of these disturbances are not the 
same for different wind-tunnels; therefore, experimental results on so-called 
natural transition obtained at different places are usually not comparable. 
We may assume that transition would take place in a similar way if levels of 
residual disturbances were smaller than certain limits. But even if this assump- 
tion is true, no one knows these limits, and it is probable that disturbances in 
existing wind-tunnels are almost always too great. 

The advantages of artificial disturbances are twofold; one is to specify the 
nature and intensity of the predominant disturbance which induces the transition, 
and the other is to improve the reproducibility of the transition process. When 
the intensity of the artificial disturbance is considerably higher than that of 
residual disturbances, the reproducibility of the transition pattern is improved 
remarkably. However, if the intensity is too high, the pattern of transition is 
completely changed from that of natural transition. The frequency and ampli- 
tude of artificial disturbances are therefore to be determined by consideration 
of these two conditions. 

A sound wave is one of the simplest and most convenient forms of artificial 
disturbance, since the frequency and intensity can be changed very easily. 
However, the excitation by sound extends uniformly over the whole flow field, 
and a point of excitation cannot be specified. Sometimes, the applicability is 
restricted by acoustical resonance of the wind-tunnel. In these respects the 


| 
ed 
= 
ut 
; 
n, 
in 
yn. 


348 Hiroshi Sato and Kyoichi Kuriki 


vibrating ribbon (Schubauer & Skramstad 1948) is a more satisfactory measure 
than sound waves. But it is difficult to apply the ribbon technique to a wake 
and, moreover, the disturbance generated by a ribbon is undesirably large in 
such highly unstable flows as wakes and jets. The transition might be induced 
by the ribbon even when it does not vibrate, and the growth of disturbances is 
completed in a small region which is not suitably large compared with the 
dimensions of the ribbon. Although the excitation by sound is not localized, 
its effect is greatest at the X-station where the spatial amplification of the 
impressed disturbance is a maximum. As a matter of fact, the velocity ampli- 
tude of the sound may be extremely small: for instance, 0-1 cm/sec for sound at 
80 db. Therefore, the induced velocity fluctuation is observable only when it is 
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Figure 28. Streamwise variation of various quantities in the transition region. 

Uy = 10-0 m/sec, Plate I. (%;)' is the maximum value at each X-station, and (%3,)* is the 

value at Y = 0. Both are plotted on an arbitrary scale. 


amplified about 100 times. Thus, it is possible to consider that the excitation is 
applied at an original point of amplification. The propagation velocity, rate of 
amplification and distribution function of spectral components in the linear 
region were determined in our experiments by the use of sound excitation. 
Wave-forms, the mean-velocity distribution and the distribution of total energy 
u? were observed without sound excitatica. 

Some of the experimental results described in the preceding sections are sum- 
marized in figure 28, which provides a general view of the mechanism of transi- 
tion. The transition region is classified into three subregions, linear, non-linear 
and three-dimensional. The linear region in which the sinusoidal velocity 
fluctuation was amplified exponentially extended from X = 20 to 40 mm, and 
this was followed by the non-linear region in which a two-dimensional double 
row of vortices was established as a stable configuration. The three-dimensional 
region started at about X = 150 mm and was succeeded by a turbulent region 
at about X = 400 mm. This result was found with Plate at U, = 10-0 m/see, 
the corresponding Reynolds number based on the length of plate being 2-1 x 10°. 
Under different experimental conditions, the relative extent of each subregion 
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was, of course, different. However, the essential features of transition were 
concluded to be identical. In order to obtain a better understanding of the 
transition phenomenon, the mechanism of transition in other flow fields is now 
to be compared with the present result. 

First, the wake of a circular cylinder is considered. Many extensive investi- 
gations have been carried out on the double row of vortices in the wake of the 
cylinder since the pioneering work by von Karman. Hot-wire measurements 
were made by Kovasznay (1949) and Roshko (1954). Roshko has pointed out 
the great difference in the flow pattern behind cylinders with different Reynolds 
number. Therefore, before we make a comparative discussion, we have to 
establish the correspondence of Reynolds number between the two types of 
wake. It seems reasonable to compare the drags of a cylinder and flat plate in 
uniform flow, which are expressed by 


in which / and d are the streamwise length of the plate and the diameter of the 
cylinder respectively. The boundary layer on the plate is assumed laminar 
everywhere. Since the drag is a measure of energy taken out of the flow, the 
velocity defect and velocity fluctuation might well be of the same level if the 
drag is equal. Thus, by putting Dyjate = Deytinaers We have 


Id = 0-38Cp 


in which R, = Ujl/v. The present experimental conditions for a plate, / = 30 em, 
U, = 10-0 m/sec, may thus correspond to d = 0-18 em and U, = 10-0 m/sec for 
a cylinder, the Reynolds number U,d/v being 1200. According to the classifica- 
tion by Roshko, this Reynolds number belongs to the ‘irregular range’ in which 
the periodic vortex shedding is accompanied by irregular velocity fluctuations. 

In the wake of a cylinder there exists no linear region in the present sense, 
that is, the sinusoidal fluctuation is of high intensity and always decays down- 
stream. The mechanism of generation of vortices is different for the two types of 
wake. Behind a cylinder there are two separated layers from the wall, which are 
laminar at first and become turbulent before the rolling-up into vortices is 
completed at Reynolds numbers below 10°. The continuity of the boundary 
layer into the wake is much smoother for the case of a flat plate. Besides these 
differences in the mechanism of generation, the double row of vortices in the 
wake of a cylinder might correspond to that in the non-linear region of the wake 
of a flat plate. Thus, in the wake-like velocity distribution, the anti-symmetrical 
two-dimensional double row of vortices is a stable configuration. The three- 
dimensional distortion of these vortices leads the gradual change into turbulence. 

The transition of a two-dimensional jet has also been found to be characterized 
by sinusoidal velocity fluctuations which grow exponentially and develop 
gradually into turbulence (Sato 1959b). Observed features of the sinusoidal 
fluctuations were in good agreement with predictions of linearized theory. 
Generally speaking, the nature of transition of a jet is very close to that of a wake 
behind a thin plate. The separated layer from a sharp edge, though it is not a 
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symmetrical flow field, exhibits a similar nature of transition to jets and wakes 
(Sato 1959a). Thus, in free boundary layers including wakes, jets and separated 
layers, the process of transition is to a large extent similar, although in the wake 
behind a cylinder there is some difference, as mentioned above. 

We now compare the transition of the boundary layer along a solid wall. The 
early work by Schubauer & Skramstad (1948) on the boundary layer along a 
flat plate has verified the linearized theory in the linear region of transition. 
Concerning the non-linear region, we can point out many essential differences 
between free and ‘fixed’ boundary layers. Modern investigations on the 
boundary layer along a flat plate indicate the importance of three-dimensional 
distortion of two-dimensional waves in the process of transition (Hama, Long & 
Hegarty 1957; Schubauer 1958; Tani 1960). The pre-existing three-dimen- 
sionality in the laminar layer—for instance, the spanwise variation of layer 
thickness—is a decisive factor in the distortion of two-dimensional waves and, 
therefore, in the onset of turbulence. Although the role of three-dimensionality 
in the transition of free boundary layers is not yet fully understood, it might be 
much less significant than in the case of fixed boundary layers. One of the most 
important features in the transition of a fixed boundary layer is the generation of 
turbulent bursts or, in other words, sudden breakdowns of laminar flow. It is 
now established that the breakdown is closely related to the complicated three- 
dimensional distortion of originally two-dimensional waves. The development of 
a small patch of turbulence has been made clear experimentally (Schubauer & 
Klebanoff 1956). These features of the non-linear region have no counterparts in 
the transition of free layers. Thus it is concluded that the linear development of 
sinusoidal fluctuations is verified in both free and fixed boundary layers, whereas 
the mechanism in the non-linear region is completely different. 

The formation of a double row of vortices in the wake is a peculiar fact. As 
mentioned in a previous paper (Sato 1959b), there might be some mechanism 
which suppresses the deformation of vortices due to a non-linear effect. The 
vortex theory, such as that originated by von Karman and modified by many other 
investigators, might explain some of the observed results. However, it is far from 
a complete explanation since it fails to describe the continuity from the linear 
region, or to predict wavelength, etc., in relation to the mean-velocity distribu- 
tion of the wake. Thus we are inclined to think that a stability calculated for 
disturbances of finite amplitude is required. It is expected that the observed 
formation of a double row of vortices is the result of a selective action which is 
closely related to the initial conditions and mean-velocity distribution. 

The three-dimensional distortion of two-dimensional vortices is initiated by 
pre-existing three-dimensionality in the flow. Recently, Benney & Lin (1960) 
have developed a non-linear theory and calculated second-order terms ¢ ssuming 
a periodic variation in the Z-direction for the primary fluctuations. From the 
present measurements of the distributions of wv? and w?in the Y-and Z-directions, 
no definite conclusion was reached concerning the validity of the theory, 
because the calculation was made for a different type of flow field and the 
assumed distribution of primary fluctuation was not realized in the present 
experiment. 
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In the wake of a thin plate, the process of three-dimensional distortion is 
gradual in contrast to the sudden breakdown in fixed boundary layers. In 
boundary layers along solid walls, there is a great difference in the equilibrium 
mean-velocity distribution of laminar and turbulent layers. This implies a more 
complicated variation of energy flow in the process of transition of fixed layers 
than that of free layers, in which thickening of the layer takes place without 
appreciable change in the shape of the mean-velocity distribution. The three- 
dimensional region of free layers is to be investigated in the future in more 
detail by determining the energy balance in the region. 


8. Conclusion 


The experimental and theoretical investigation of the wake developed behind a 
flat plate placed parallel to the uniform flow indicated the following conclusions. 

1. The observed mean-velocity distribution in the laminar region of a wake 
is in coincidence with Goldstein’s calculation. 

2. At acertain Reynolds number, the transition region is classified into three 
subregions: (1) the linear region, (2) the non-linear region, and (3) the three- 
dimensional region. 

3. The frequency of sinusoidal velocity fluctuations found in natural transition 
is proportional to the 3 power of free-stream velocity. The fluctuations corre- 
spond to the disturbance with maximum amplification according to the linearized 
stability theory. 

4. The sinusoidal velocity fluctuations induced by artificial excitation are 
amplified exponentially in the linear region when the frequency of excitation 
lies inside the unstable zone of the stability diagram. The observed propagation 
velocity, rate of amplification and amplitude function of fluctuations at various 
frequencies are in good agreement with the predictions of linearized theory. 

5. Vortex models, a single row of vortices in the linear region and an anti- 
symmetrical double row in the non-linear region, are useful for explaining many 
of the experimental results. The double row of vortices is similar to that found 
in the wake of a circular cylinder. 

6. In the non-linear region, an increase of intensity of the artificial excitation 
results not only in the increase of induced velocity fluctuation of identical 
frequency, but also in a decrease in the amplitude of the natural velocity 
fluctuations. This ‘suppression effect’ is remarkable when the frequency of 
the artificial disturbance is close to that of natural fluctuation. 

7. In the three-dimensional region, the double row of vortices, which is two- 
dimensiona! in the non-linear region, is subject to three-dimensional distortion. 
The velocity fluctuation gradually develops into turbulence without being 
accompanied by sudden breakdown or turbulent bursts. 


The authors appreciate the stimulating discussions occurring throughout the 
whole course of the investigation with members of the Boundary Layer Research 
Group in Japan, which is directed by Professor Itiro Tani. Thanks are extended 
to Mr Y. Onda and Mr O. Okubo who helped to carry out the experimental work. 
The research was financially supported by the Ministry of Education of Japan. 
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The behaviour of a natural convection plume above a line fire is studied both 
theoretically and experimentally. In the theoretical treatment, a turbulent 
plume above a steady two-dimensional finite source of heated fluid in a uniform 
ambient fluid is investigated. By the use of the lateral entrainment assumption, 
a quadrature solution has been obtained for each of two separate ranges of a 
source Froude number, / > 1 or F < 1. In neither of these cases can the finite 
width line source be accurately represented by an equivalent mathematical 
line source at a lower level. Only the special case, F = 1, can be so represented 
and its solution is discussed. 

In the experimental treatment, hot gases, resulting from the burning of a 
liquid fuel in a long channel burner, are driven upwards by buoyancy and gradu- 
ally cooled down by the entrainment of ambient air. The average temperature 
along lines parallel to the channel burner was measured by a piece of resistance 
wire. For the case of a non-luminous diffusion flame, the effective radiation loss 
to the surroundings was assumed to be negligible, and, by a comparison of the 
energy flux supplied from the fuel and the energy flux contained in the plume, 
the characteristic turbulence entrainment coefficient is determined. By the alter- 
nate use of either an absorbing or a reflecting surface for the table-top surround- 
ing a luminous flame, a measurement was made of the energy radiated from the 
flame that was intercepted by the fire surroundings and subsequently returned 
to the buoyancy plume by heating the ingested air. These measurements agree 
with estimates computed from such data as are available. The experimental 
results relating to the behaviour of the convection plume agree closely with the 
theoretical predictions in all cases. 


1. Introduction 


A line fire exhibits a complex interaction of diverse phenomena. Hot combus- 
tion products are driven upwards by natural convection, thus inducting air and 
fuel for further combustion. Heat is transferred downward by various mechanisms 
to maintain the combustion temperature in the fresh fuel and air and, if the fuel 
is liquid, to continue the evaporation process. Simultaneously, interchange of 
energy between the fire and the various elements of its environment also occur. 
One of the effects of this phenomenon is the instability of the inflowing air, 
which results from heat transferred from luminous flames to the surrour.Jing 
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surface by radiation and thence to the inducted air by conduction and convection 
as noted by Emmons (1959). 

Vo attempt is made to solve the flame-fuel evaporation problem, but rather 
attention is given only to the natural convection above the flames. Only the 
turbulent flow régime is considered. Experimentally, measurements are made 
above the level of the flames while theoretically the whole fire process is replaced 
by a horizontal source of heat, momentum and energy, infinitely long and of 
finite width. In the theory presented these may be varied independently. Since 
the experiments were conducted in the laboratory, the plumes were only a few 
metres high and atmospheric variations were negligible. The present work, there- 
fore, extends previous work by the inclusion of arbitrary fire width, and arbitrary 
momentum and energy sources, rather than the line fire only. 


2. Theoretical investigation 


It is assumed that turbulent flow is fully developed and thus, as a consequence, 
molecular transier mechanisms are negligible relative to the turbulent transfer. 
It is also assumed thai local density variations are everywhere small compared 
to some reference density in the field. Although natural convection of this 
nature is produced by a source of heat, it is the buoyancy rather than the thermal 
properties of the flow which is fundamental to the phenomenon. Therefore, 
although the buoyant force due to density difference is sufficiently great to 
produce vertical acceleration, the corresponding variation in the mass density 
of the fluid undergoing acceleration is sufficiently small, in comparison with the 
density itself, to be neglected in the continuity and energy equations and the 
inertia term of the momentum equation. 

In addition, we make the ‘boundary-layer’ assumption that transverse 
accelerations are small in comparison with those in the vertical direction, and 
that turbulent mixing in the vertical direction is small in comparison with that 
in the horizontal direction. A consequence of this is that the pressure has essen- 
tially no horizontal variation. 

Let the origin O represent the location of the centre of a two-dimensional 
finite source of width b) (accurately defined by equation (8) applied at 2 = 0). 
The local mean component of velocity in the vertical direction x is denoted by w, 
and that in the y-direction. (in the horizontal plane, normal to the fire) is 
denoted by v. 

Under the foregoing assumptions the fundamental equations of motion reduce 
to the equation of continuity 


Cu Cv 
= 0, (1) 
Cu Cy 


the simplified equation of conservation of vertical momentum 


Ou dw Ay 


(2) 
py py Cy 
and the simplified equation of conservation of energy 
o(AT) | c(AT) 1 ow (3) 
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Here p, is the mass density of undisturbed ambient fluid, and 
Ay = gAp = 9(P1.—P) 
is the local buoyancy, where p is the local mass density, y the local weight den- 
sitv, and g the gravitational acceleration. Also, T= —p,uyv, is the Reynolds 
stress, u;, v, being fluctuating velocities in the x- and y-directions respectively; 
AT = T —T,, where T is the local time-mean temperature, and 7, the tempera- 
ture of undisturbed ambient fluid. Lastly, C,, is the specific heat at constant 
pressure, and w= —p,C,,v,AT; is the eddy heat transfer, where AT; is the 
fluctuating temperature difference. 
Since v = 0 at y = 0 by symmetry, integration of equation (1) with respect to 
y gives 


d 
A udy +v(x, = 0. (4) 


This equation states that the increase of ascending mass in a heated plume is due 
to the lateral entrainment of the ambient fluid. 

Since no work is applied at either y = 0 or y = © and the flow is symmetrical 
with respect to y = 0, therefore 7 = 0 at both y = 0 and y = o. Integration of 
equation (2) with respect to y by parts, with use of equation (1), gives 


d [2 Ay 
0 Pi ©) 


This equation states that the increase of vertical momentum of ascending mass 
of heated fluid is due to the buoyancy effect caused by the density difference. 
Since no heat is added at either y = 0 or y = © and the flow is symmetrical 
with respect to y = 0, therefore w = 0 at both y = 0 and y = oo. By the assump- 
tion that the ambient fluid is undisturbed, A7' = 0 at y = «. Integration of 
equation (3) with respect to y, using equation (1) and the assumption of small 


density variation, gives d f= 
al uAydy = 0. (6) 
dx 0 


This equation states that the buoyancy flux of an ascending mass of heated fluid 
is conserved. 

When a stream of fluid is in contact with another stream, the eddies which 
cause transfer of matter between them are characterized by velocities propor- 
tional to the relative velocity of the two streams. This can be seen from dimen- 
sional considerations, since if the mutual entrainment is turbulent the only 
quantity determining the motion is the relative velocity of the two streams. This 
was also shown by the experiment of Kuethe (1935) on the mixing at the edge 
of a jet. The rate at which the bounding edge of a heated plume expands and 
absorbs the surrounding fluid into the plume may similarly be assumed propor- 
tional to the velocity at that level. As first suggested by Taylor (1945) and later 
adopted by Morton, Taylor & Turner (1956) and Morton (1959), the rate at 
which fluid is entrained into the heated plume is taken as proportional to the 
vertical velocity on the axis of the plume, thus 


v(x, 00) = au(x, 0), (7) 


where « is the entrainment constant. 
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Rouse, Yih & Humphreys (1952) have obtained velocity and temperature 
measurements of heated plumes above a line of small gas flames designed to 
simulate a line heat source. Their results show the Gaussian profiles 

exp (—32y?/7") and exp(—41y?/2") 

for the velocity and temperature measurements, respectively. Since the entrain- 
ment constant & is by definition associated with the ve!:ty profile rather than 
the temperature profile, it will be appropriate to seek: » sunilarity solution for 
which x measures the rate of flow into a heated plume, with velocity profile 
characterized by a horizontal length-scale b, and with an associated buoyancy 
profile of the same shape but with a length-scale Ab. « and A are assumed to be 
universal constants. 

Consider a two-dimensional heated plume generated from a finite source in 
an incompressible environment. The plume will be assumed to have Gaussian 
profiles of time-mean vertical velocity and time-mean buoyancy: thus 


u(x, y) = u(x) exp (—y2/b%), (8) 
where w(x) = u(x, 0) is the time-mean vertical velocity on the axis of symmetry 
of the plume, and Ay(x,y) = Ay(x) exp (— y2/A2b2), (9) 


where Ay(v) = Ay(v, 0) is the time-mean buoyancy on the plume axis. The finite 
source will be characterized by = bo, u(O) and Ay(0) = Sub- 
stituting equations (8) and (9) into equations (4), (5) and (6), we have the equa- 
tions of continuity, conservation of vertical momentum and conservation of 


energy 
= 2) (11) 
V1 
d 
0, (12) 
subject to the boundary conditions 
b=b, u=uU, Ay=Ay, at x=0. (13) 
With the following transformations designed to remove all unnecessary 
coefficients, x’ = (2/3) xa/bo, 
b’ = b/bo, 
= (14) 


we have, from equations (10) to (13), the equations 


21al 


= = 2", (15) 
d 
Lu = (16) 
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subject to the boundary conditions 


w=F, p=l at =90, (18) 
, (2%flal ) 


We note that this group can be written . 


7 A Ayo, (gbo)? 
a source Froude number of the usual form w/(gb,)* moditied by factors (2/7)! 
to account for the Gaussian distribution, («/A)! to account for the difference 
between velocity and buoyancy profiles, and (Ayp/77,)? to account for the actual 


effective density difference. This Froude number plays an essential role in what 
follows. Equation (17) can be immediately integrated to give 


u'b'p’ = F, (19) 
by the use of which equation (16) becomes 


d F ) 


To solve equations (15) and (20) introduce the ‘mass flux’ and ‘momentum 
flux’ as two new variables 


(21) 
dN M 
Thus we get (22) 
dM FN 
dx’ (23) 


subject to the boundary conditions 
M=F* = 6. (24) 
By eliminating 2’ and integrating, we find that 
FN? = +c, (25) 


where c = F4(1— F?). We note that both mass and momentum flux are positive 
and by equations (22) and (23) increase with increasing x’, although ¢ may be 
zero or have either sign. 

The position in the plume is now found by integrating equation (22). i.e. 


CN AT y 
v=| NdNn | NdN (26) 


r (FN3—c)s 
It is convenient to consider separately three cases, c < 0,c = 0,c > 0, and 
arrange the integral for graphical or numerical evaluation. The three ranges 


correspond to F < 1, F = l,andF > 1, respectively. 
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Case I. F = 1, the line source 


For F = 1, ¢ = 0, equations (25) and (26) give 


and N =f, (28) 
which together with equations (19) and (21) give 
1 
w=1, =2'4+1, p’ 29 
(29) 
These equations are plotted in figure 1. 
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Dimensionless vertical distance 2’ 
Ficure 1. Graphical representation of the theoretical solution for a plume issuing 


from a neutral source, F = 1. 


If we introduce the conserved buoyancy flux across any horizontal cross- 


section 
—o 1 


we have, from equation (29), the familiar forms 


b= ax+bo, (31) 
v1 
Ay (1+A2)3 Qi 2a 


If some arbitrary point in the plume were regarded as the source at which the 
convection started, we would compute a source Froude number of F = 1 inde- 
pendent of location. Thus, the above solution describes the convection above a 
line source 6, = 0 in which case the velocity uw = const., the plume ‘width’ b 
is proportional to x, while the buoyancy Ay/y, is inversely proportional to 2. 
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Furthermore, for a finite width line source, the whole convection column can 
be regarded as arising from a line source situated a distance a) = 7b)/2x below 
the real source. We note that the equivalent line source has a finite buoyancy 
flux if by Ayo is constant; i.e. if the density difference becomes infinitely great as 
the source width goes to zero. 


Case II. Plume issuing from a restrained source (F < 1) 
The integral of equation (26) is expressible as an incomplete beta function. 
However, direct numerical integration is more convenient. For this purpose the 
parameters are removed by setting 


NY = fy, (32) 
with p= () = F(1—F?)}. 
There results 

v sdp 

(33) 
9\—1 N 

where y = (1—F?)-3 7 | (34) 


= 


These integration limits are both greater than 1, the upper limit being larger 
than the lower limit. Since for the Froude number limits of this case F = 0, 
F = 1, v has limits of 1 and oo, we define 


vdv 
10) = (35) 


and the height scale in the convection plume is given by 
= {F%(1—F*)}} (36) 


Thus the actual plume may be regarded as a section of the plume arising from 
a virtual source at a distance {F2(1—F?)}} J,(v)) below the real source. The 
strength of the real source is specified by its ‘mass flux’ Ny = F, while the ‘mass 
flux’ increases along the plume as N = F(1— F?)! p increases. 

From equations (19), (21), (25) and (32), the physically interesting quantities 
are given as follows: 


plume half-width = {F?(1— 
p 


The integral J,(v) was evaluated graphically and the physical plume charac- 
teristics are plotted in figure 2. The real plume may be visualized from this figure 
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by imagining the actual source located at a distance to the right of the origin 
given by 
0 

38 
while other higher plume points appear to the right of the actual source. The 
plume shape is particularly interesting. For areal source with F < 2- the plume 
width first decreases and then increases again. Furthermore, any plume at ¢ "eat 
height behaves as though it arose from a line source located at a distance (in 
the scale of equation (37)) of 0-7 below the virtual origin. 


a 


Ww 


1/(1— p’ 


Dimensionless vertical velocity u’/F} 


vy 


Dimensionless half-width b’/[F2(1 — F?)]} 
Dimensionless reciprocal of buoyancy 


2 
Vv Fi 
=a 
0 1 2 3 4 


Modified distance above virtual source J,(v) 
Ficure 2. Graphical representation of the theoretical solution for a plume issuing 
from a restrained source (see equation (33)), F < 1. 
Case III. Plume issuing from an impelled source (F > 1) 


We again remove the parameters from the integral of equation (26) by setting 


N = pf, (39) 
where now = F(F?-1)3. 
The resulting integral for graphical evaluation is 
vdv 
a’ = {F(F2-1 40 
(40) 
with v= (41) 
Vo = (F?- 1)+4, 


where 0 < Vy < V < ©. 
Fer this case we use zero as the lower limit and define 


(42) 
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rin from which the position in the plume is computed as 
8) a’ = {F*(F2—1)} {1,(v) — 1,(v9)}.- (43) 
By equations (19), (21), (25) and (39), the physically interesting quantities are 
he obtained: 
plume half-width 6’ = {F?(F?—1)}» ———,, 
ne (v3 + 1)3 j 
at 3 
buoyancy 


1/(F2—1)8 p’ 
Dimensionless vertical velocity u’/F4 


vy? b 
+1) 
y+ 


Dimensionless half-width 
Dimensionless reciprocal of buoyancy 


0 1 2 3 4 


Modified distance above virtual source J,(v) 


FicureE 3. Graphical representation of the theoretical solution for a plume issuing 
from an impelled source (see equation (43)), F > 1. 


Curves of these plume properties are given in figure 3. Again the origin of the 
distance scale may be regarded as a virtual source with the real source located 


at a distance 
—1)h, 45 
(45) 
and the real plume given by the curves from the real source to infinity. The plume 
grows in size rapidly at first and then more and more slowly. In this case too, the 


) plume at great heights becomes like the plume from a line source located at a 
distance of 0-7 below the virtual source. 
) We note that both cases II and III approach the case of F = 1 as height 


increases. F = 1 corresponds to a balance between inducted mass increase and 
buoyancy which just maintains the vertical velocity constant. If the plume 
velocity is relatively too small, i.e. F < 1, the plume grows slowly or even con- 
tracts to raise the local Froude number. If, on the other hand, the plume velocity 
is too high, as in a heated jet, i.e. / > 1, the plume grows more rapidly, decreasing 
the velocity. 
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3. Experimental investigation 


The channel burner, 0-564in. wide, (-282in. deep and 78in. long, is cooled 
on the sides by a water-jacket, which controls the fuel consumption rate, and 
is connected by inlets underneath to a metered steady fuel supply system. On 
either side of the channel burner there is an 18in. wide asbestos side plate set 
flush with the sides of the channel burner to prevent any entrainment of air from 
below. Against each end of the channel burner stands a vertical asbestos end 
plate to prevent the entrainment from the ends. A view of the apparatus is 
shown as figure 4 (plate 1). Except for cases with very low fuel consumption rates, 
the flame from the channel burner fluctuates both with time and location along 
the length of the channel burner. The average temperature at different heights 
was measured with a 40in. piece of pure nickel resistance wire 0-002 in. in dia- 
meter suspended parallel to the channel burner. Two fuels were chosen, acetone 
and methyl alcohol, because they differ considerably from one another in radia- 
tion from the flames. Two different kinds of surface were chosen for the side 
plates, aluminium foil surface and bare asbestos surface, because they react 
very differently to radiation. 


(A) Analyses of asymptotic behaviour of plumes 


For either methyl alcohol or acetone, Gaussian profiles of the natural convection 
plume are found to prevail at all measurable heights and the corresponding values 
of the characteristic buoyancy half-width are thus computed. For each plume 
the values of the slopes d(y,/Ay)/dx and d(Ab)/dx were obtained graphically from 
the measured centre-line buoyancy and Gaussian width respectively and are 
found to approach their respective asymptotic values. 

These asymptotic values are given in a convenient form by equation (31) 


dx lim Qi | 


| 
dx lim 7 : 


The energy flux Q in this equation could be computed from the rate of fuel use 
and the fuel properties if it were not for radiation losses. If radiation is negligible. 
the energy flux from the fuel Q is 


(46) 


= 47 
Q (47) 


where q¢ represents the mass consumption rate of fuel per unit length of burner, 
H the heat of combustion per unit mass of fuel, and dh the heat removed by the 
cooling water per unit mass of fuel burned. 


(a) Convection plumes above a non-luminous diffusion flame of methyl alcohol 


Since radiation is relatively unimportant in the process of heat transfer from 
a non-luminous flame, the surface condition of the side plates would be expected 
to have relatively little effect on the behaviour of the convection plume above 
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the flame. Experiments with asbestos board side plates and bright aluminium 
foil side plates gave nearly identical convection plumes with the methyl] alcohol 
flames. Thus the heat released from chemical reaction minus the amount removed 
by cooling water is carried away by the convection plume and Q equals Q. 
Therefore the values of the assumed universal constants « and A can be computed 
from equation (46) and the asymptotic data. We find « = 0-16, A = 0-9, which 
check very well with values computed from results of velocity and temperature 
measurements by Rouse ef al. (1952) from natural convection plumes above a 
simulated line buoyancy source. Furthermore, the value of « also checks very 
well with a value computed from velocity measurements of a turbulent two- 


dimensional free jet by Reichardt (1942). 


(b) Convection plumes above a luminous diffusion flame of acetone 

Since in a luminous flame carbon particles radiate energy to the surroundings, 
the value of Q in the convection plume is usually smaller than that of Q. The 
difference Q—Q is the net radiation loss from the flame-plume system to the 
surroundings. The results of experimental measurements with asbestos and with 
aluminium side plates are shown in figure 5. 


Fractional radiation loss (Q —Q)/Q 


retic edicti = 
Acetone based on a simplified radiation model ‘O@ 8 = 
0-40- 40-40 op 
23 
Experimental results of iit 
zee 


Ratio of radiation losses with absorbing 


0 0-002 0-004 5006 0008 
Fuel consumption rate (g see~! 
FicureE 5. Radiation loss computed from measurements of the plume above a luminous 
diffusion flame of acetone: [], with bare asbestos side plates; ©, with aluminium foil- 
covered side plates (estimates based on average temperature of carbon particles and 
emissivity measured by Tankin 1960). 


It is readily seen that the net loss to the surroundings is lower for the case with 
bare asbestos side plates, than for those covered with aluminium foil. This can 
be explained by the difference in the way radiant energy is dealt with by the 
surface of the side plates in the two cases. The bare asbestos surface has a very 
high emissivity (¢ = 0-96) and thus absorbs almost all of the radiant energy from 
the flame intercepted by it. The asbestos surface thus heated will give up heat by 
conduction to the adjacent air, which is then convected into the flame-plume. 
On the other hand, the aluminium foil surface has a very low emissivity 
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(¢ = 0-087) and will reflect almost all of the radiant energy from the flame 
intercepted by it. Thus all radiation from the flame not directly intercepted by 
the fuel is permanently lost from the convection column. 

To estimate the flame radiation loss, use is made of measurements by 'Tankin 
(1960) with a bolometer of the average temperature of flame carbon particles, 
and of the emissivity of a luminous diffusion flame of acetone. A simplified line 
radiation source was used to estimate the radiation loss of a flame of approxi- 
mately the same brightness, over the channel burner. The results, shown as solid 
points in figure 5, are in good agreement with findings of the plume measurements. 
The ratio of radiation losses with absorbing and with reflecting side plates based 
on the same simplified radiation model, is also shown in figure 5 as the broken 
curve which checks quite well with findings of the plume measurements at all 
values of fuel consumption rate. 


(B) Comparison of theoretical and experimental convection plume results 


It is predicted by the theoretical analysis that the behaviour of a plume above a 
finite source depends on a parameter F which is a function of the source charac- 
teristics, by, uy and Ay,/y,. It is then desirable to determine the values of F’ by 
comparing the theoretical predictions with the experimental data. 

It is seen from figures 2 and 3 that the solutions for the variation of plume half- 
width with vertical distance and the location of the actual source above a virtual 
source, when all dimensionless quantities are converted back to physical quanti- 
ties, depend only on two of the three independent characteristics of the source, F 
and 6). For an origin assumed to be at the mean flame height, experimental 
data on the plume half-width are plotted against theoretical predictions based on 
different pairs of trial values of F andby. By the principle of least squares, the best- 
fit values of F' and by are determined for each case. The corresponding best-fit 
value of Ay,/y, is then determined by comparing the theoretical predictions 
and experimental data on the variation of buoyancy with vertical distance. The 
values of F thus found for an origin assumed to be at the mean flame height 
are listed in table 1. 

The plumes from a methyl alcohol flame are seen to lie in the region of F < 1, 
while the plumes from an acetone flame are found to lie in the region of F < 1 for 
low burning rates and F > 1 for high burning rates. Very good agreement has 
been found between theoretical predictions and experimental measurements for 
ail the experiments performed as shown by the sample comparisons made in 
figures 6 and 7. 


(C) Interpretation of experimental results with the introduction 
of a simplified flame model 


The rigorous analysis of a flame is extremely complicated. However, if the 
main interest is in the study of the convection plume then the flame itself serves 
merely as a supplier of heated fluid to the source at the mean flame height for the 
convection problem. Therefore, if we introduce a simplified flame model with 
uniform temperature and vertical velocity distributions along a mean flame 
contour as shown in figure 8, the characteristics of the source can be determined. 
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Based on such a simplified flame model, the energy flux from the flame is 
V1 


where 6, is the half-channel width of the inside burner. The energy flux contained 
in the plume, from equation (30), is - 


Q = (at (49) 

Fuel Side plates q (g see cm) F 
Acetone Bare asbestos 0-00110 0-79 
0-00163 0-98 
0-00191 1-33 
0-00268 
0-0055 3°56 
0-0093 5°85 
0-0138 8-67 

0-0170 10-7 
Acetone Aluminium foil 0-00100 0-65 
0-00134 0-85 
0-00157 0-98 
0-00254 1-49 
0-0041 2-14 
0-0070 3-74 
0-0103 5:37 
0-0127 6-75 
0-0165 8-82 

0-0187 10-3 
Methyl Bare asbestos 0-000716 0-39 
alcohol 0-000885 0-48 
0-00158 0-85 
Methyl alcohol Aluminium foil 0-0074 0-40 
0-:00158 0-85 


TasBLe 1. Values of Froude number at the effective source obtained from 
plume measurements. 


Comparing equations (48) and (49), we have 


2 
(50) 
Eliminating uw) and b) among equations (18), (49) and (50), we then have 
1 ( V1 ) 
51 


The values of F of the source computed from equation (51) with such data as 
available check closely with those determined from plume measurements as 
shown by the sample comparison in figure 9. 
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4. Conclusions (a) Theoretical conclusions 


A quadrature solution for a two-dimensional turbulent natural convection plume 
above a finite source has been obtained with the assumptions of lateral entrain- 
ment and similar Gaussian velocity and temperature profiles. A characteristic 
parameter F’ based on the characteristics of the finite source, has been found to 
describe the behaviour of the plume above that source. Plumes for cases where 
F < land F > 1 behave differently but both approach asymptotically the case 
where Ff’ = 1. The solution for the case where F = 1 is identified as the solution for 
a plume above a line buoyancy source situated at a lower level. 
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FicurE 6. Example of comparison between theoretical and experimental results for 
the buoyancy in the plume above an acetone flame, Ff > 1. 


q q 
(g (g sec“! 

F 
fr © 0:00191 1-33 + 0-00254 1-49 
© 000268 1-83 Side plates 00041 2-14 
Bare asbestos © 0:0055 3-56 covered + 0-0070 3-74 
side plates ©- 0-0093 5:85 with 00-0103 5:37 
© 00-0138 8:67 aluminium =0-0127 6°75 
0-0170 10:70 foil 60-0165 8-82 
0-0187 10-33 


(b) Experimental conclusions 


(i) A Gaussian profile was found to prevail at all measurable heights. 

(ii) Values of the entrainment constant « and the auxiliary constant A obtained 
from methyl alcohol burnings agree closely with those evaluated from the pre- 
viously established data. 
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| FigurE 7. Example of comparison between theoretical and experimental results for the 
half-width of the plume above an acetone flame, F > 1. The theoretical curve is drawn 
in a continuous line. 
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for the convection plume. 
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(iii) The radiant heat loss from a luminous flame to its surroundings deter- 
mined from measurements on convection plumes agrees with estimates computed 
from such data as are available. The effect the surface condition of the side plates 
has on the behaviour of the plume from plume measurements checks very well 
with that from an estimate based on a simplified radiation source. 

(iv) Experimental results for the natural convection plumes above a line 
fire are found in all cases in good agreement with theoretical predictions for the 
two-dimensional natural convection above a steady finite source of heated fluid. 
An interpretation is provided by the assumption of an effective source based on 
a simplified flame model. 


8k 


Froude number F 


0-002 0004 0006 0:008 0010 0012 0014 0016 0018 0020 

Fuel consumption rate ¢g(g sec-! 


Ficure 9. Example of comparison between the computed Froude number based on a 
simplified flame model and the experimentally determined Froude number of an effective 
source for the convection plume above an acetone flame: @, with bare asbestos side plates; 
©, with aluminium foil covered side plates. The two theoretical curves are drawn as 
continuous lines, the respective curve being the one closest to the experimental results. 
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Experiments on the penetration of an interface 
by buoyant thermals 


By J. M. RICHARDS 


Department of Meteorology, Imperial College, London* 
(Received 3 January 1961 and in revised form 11 April 1961) 


Isolated masses of a dense aqueous solution (i.e. thermals) were released at the 
surface of a freshwater layer overlying one of salt water. While the whole of a 
thermal remained in the upper layer, the equations z = n,r, 2? = k,t, with 
k, = C,n}(Mg|p)}, were obeyed, where z is the distance travelled, 27 the width of 
the thermal, ¢ the elapsed time, 17 the mass excess (the difference between the 
masses contained within and displaced by the thermal while in the upper layer), 
and p the density of the displaced fluid. », was constant for any one thermal, 
but varied between thermals over the range 1-9 to 7-5. C, had roughly the same 
value (0:73) for all thermals. 

The behaviour after the leading extremity of the thermal (the ‘front’) entered 
the lower layer depended on the value of the parameter S = VAp/M, where Ap is 
the density difference between the upper and lower layers and V is the volume of 
the thermal when the widest part is at the level of the discontinuity. It was found 
that Y=0 if S>/ (‘weak’ thermals), and Y =0-95-48 if 01<S<f 
(‘strong’ thermals), where Y is the fraction of the mass of the substance released 
which penetrated indefinitely into the lower layer. The constant / was approxi- 
mately equal to 1-90. 

In weak thermals, the equation z—s = a,(t—t,)? was obeyed while the front was 
in the lower layer, until the culminating point z = s was reached at t = ¢,. The 
acceleration a, was always negative, and constant for any one thermal, but varied 
between thermals. Also for weak thermals, 2 = C, V!/S, where z is the distance 
from the interface to the culminating point and C, is a constant. C, was approxi- 
mately equal to 3-5. For strong thermals, the distance travelled while the front 
was in the lower layer obeyed the equation z? = k,t. The origins of z and ¢ 
usually differed from those found in the upper layer, and generally k, + k,. 


1. Introduction 


In recent years, early attempts to describe the convection processes in the 
atmosphere which are responsible for the initiation and growth of cumulus clouds 
have led to ‘bubble’ and later ‘thermal’ models of buoyant convection (Bjerknes 
1938; Stommel 1947; Davies & Taylor 1950; Scorer & Ludlam 1953; Scorer 1958). 

In the experiments on the thermal model described by Scorer (1957), masses of 
marked salt solution were observed as they moved from rest at a free surface, 
under the action of buoyancy forces alone, into initially still water of uniform 

* Now at Department of Electrical Engineering, Loughborough College of Technology. 
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density. It was seen that most of the marked fluid entered a well-defined region, 
or ‘thermal’, which moved through the surroundings in the direction of the 
buoyancy force. After an initial acceleration the shape of the thermal remained 
nearly constant, but its size increased continually and the rate of growth and 
speed of advance decreased continually. 

In the experiments reported here, masses of opaque solutions were released at 
the surface of a comparatively shallow layer of water superimposed on a deep layer 
of sodium chloride solution of uniform density. As the lowest part of the thermal 
(‘the front’) approached the interface, the varying pressure distribution over the 
interface created waves, but the maximum displacement of the interface was 
always very small compared with any linear dimension of the thermal cloud while 
this remained entirely within the upper layer. 


2. Motion in the upper layer 
Dimensional analysis 


The greatest linear dimension of every thermal was initally much less than the 
depth of the water layer, and we may therefore assume that the effect of the 
interface upon the early part of the motion was negligible. The accelerated part 
of the motion was brief, and will not be discussed. It will be assumed that the 
remainder of the motion may be sufficiently described by the following physical 
variables: the speed of advance of the front, dz/dt; the rate of increase of the 
maximum horizontal dimension of the thermal, 2dr/dt; the total buoyancy force, 
Mq; the density of the surrounding fluid, p; andr. As usual, it is assumed that the 
fluid density may be taken as substantially equal to p everywhere except for the 
evaluation of buoyancy forces, and that the fluid viscosity is irrelevant to the 
motion. 
A straightforward dimensional analysis then gives 


r(dr/dt) = (Mg/p)}, (1) 
where n, is defined by dz/dt = n,(dr/dt). (2) 


It was observed (see below) that n, was approximately constant only within 
experiments. Since /qg/p was also constant in each experiment, (1) indicates that 
r(dr/dt) was constant within experiments. The Reynolds number r(dr/dt)/v was 
therefore constant within experiments; and if the distribution of fluid velocities 
in each experiment is assumed to have been similar at all times, it follows that the 
circulation K round a circuit passing along the axis of an axisymmetric thermal 
and returning in the surrounding fluid, having the same dimensions as r(dr/dt), 
was also constant within experiments. The dimensionless parameter A2p/Mg was 
therefore constant within experiments, but may be supposed to have varied 
between experiments as a result of variations in the initial conditions. 


Vortex-ring model 


A kinematic model of thermal motion has been advanced by Turner (1960). 
According to this model, a thermal may be regarded as a buoyant vortex ring 
having axisymmetric distributions of velocity and buoyancy. It is assumed that 
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these distributions are similar at all stages of the motion and that the circulation 
K, already defined, is constant within thermals. Turner showed that, for such 


a ring, x, = (K°p/Mg) mR, (3) 


where £ is the radius of the circular axis of the vortex system (Lamb 1932, § 162), 
x, is the distance of the plane of the circular axis from a fixed virtual origin, and 
m is a dimensionless function of the vorticity distribution. Since, by the 
similarity assumption, R=ar, z= ba, (4) 


where a and b are constant within thermals, (3) may be written 

z = abm(K*p/Mq)r; 
that is, n, = abm(K?p/Mq). (5) 
If the term abm is constant between thermals, this equation shows that n, is 
directly proportional to (K*p/Mq). 


Experiments 


Thermals having various mass excesses were released from a hemispherical cup. 
The cup was pivoted about a diametral axis which lay slightly above the free 
surface of the water in an experimental tank. This surface was initially coplanar 
with the surface of the fluid contained in the cup. Release was effected by rapidly 
inverting the cup about its fixed pivots. 

The mass excess in the cup was found by measuring the volumes and densities 
of the solutions used. Standard gravimetric draining-error determinations showed 
that the mass excess released was always approximately 99 °% of the mass excess 
in the cup. The corresponding values of (z — z)), 27 and (t —t,) were measured from 
enlarged images of cine photographs (2, is the vertical displacement of a fixed 
mark from the virtual origin z = r = 0, and t = f, at an arbitrary fixed instant). 
The separation of the images of two fixed marks at a known distance apart was 
measured to determine the scale of the photographic image. A graph of the 
observed values of (z —z)) was plotted against 27, and in each case a straight line 
fairly represented all except the earliest observations. The corresponding value of 
n, was found as twice the slope of this line (see (2)), and z) was determined as 
minus the intercept of the line on 2r = 0. 

From (1) and (2), we have 

z2=k,t, where hk, =f(n,)(Mg/p)}. (6) 
Graphs of z? against (¢ — t)) were plotted from the observations of each experiment. 
Straight lines fairly represented all except the earliest observations. This con- 
firmed that the interface had little effect on the motion in the upper layer. The 
value of f(n,) for each experiment was found from the slope of the corresponding 
graph and the measured value of J, use being made of (6). The pair of graphs 
relating to one particular experiment are shown in figure 1. 

Figure 2 shows a graph of the corresponding values of f(n,) and n, for each 
experiment. It will be seen that a wide range of values of n, was obtained, and 
that an equation of the form 

f(m) = Cyn}, with CO, = 0-73, (7) 
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Ficure 1. An example of motion in the upper layer. 
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FIGURE 2. Relation between the angle of spread and rate of progress of thermals moving 
from a free surface into unstratified surroundings. The curve shown by the broken line is 
that suggested by Scorer (1957). ©, Richards; x, Scorer. 
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fits the observations reasonably well. The value of C, was found as the reciprocal 
of the slope of the line in figure 2. Figure 2 also shows the observations reported 
by Scorer (1957). These are consistent with (7) and the present observations, but 
the latter are inconsistent with the empirical relation, f(n,) = 2-9n} (shown in 
figure 2 by a broken line) which was suggested in Scorer’s paper. 


Ap 
M (g/c.c. 4 ky h 
Expt. (g) x 10-8) (1.) ¥ (em?/sec) ny c (cm) 

7 18-5 5-74 1-78 0-49 501 2:8 0-76 24 

8 11-2 5-74 3-40 0-09 501 2-9 0-45 29 

9 14-9 1-27 6-98 0-70 362 2-0 0-79 24 
10 22:3 1-27 3°86 0-84 0-69 
11 15-0 1-27 3-41 0-83 829 4-7 0-85 39 
12 45:3 1-27 4-20 0-89 658 2:7 0-75 26 
13 15-0 1-27 2-84 0-77 421 2-4 0-60 27 
14 15-0 5-80 1-20 0-74 1085 5-6 0-60 41 
15 11-3 5-80 1-19 _- 714 4:8 0-84 37 
16 15-2 5-80 2-68 0-77 912 4-9 0-78 41 
21 15-1 3-92 3-03 0-60 496 2-6 0-72 27 
22 15-1 3°84 4-23 0-38 756 3-9 0-87 35 
23 15-1 3-65 1-66 0-67 1080 5-4 0-66 40 
24 15-1 3-68 3-91 0-42 565 3-4 0-75 33 
25 15-2 3°95 3°23 0-51 511 2-7 0-90 26 
26 11-3 3-74 2-68 0-49 468 3-1 0-85 29 
27 9-47 2-10 1-89 0-69 652 4-0 0-79 32 
28 4:83 0-81 1-99 0-71 459 3:8 0-76 31 
29 4:87 1-42 2-07 0-48 475 4-5 0-60 38 
30 6-04 1-07 2-60 0-73 326 3-1 0-76 26 
31 4-83 1-51 2-03 0-64 319 3:7 0-84 30 
32 15:1 1-51 3-46 0-77 464 2-9 0-66 25 
34 12-0 4-93 2-05 0-53 891 4-6 0-57 37 
35 10-6 5:83 2-98 0-07 551 3:8 0-89 32 
36 10-6 5:97 2-99 0-21 747 4:9 0-74 39 
38 9-01 5-44 3-16 0-04 721 4-4 0-77 41 
39 9-01 4:27 1-33 0-57 1048 7:3 0-63 45 


TABLE 1. Measurements on strong thermals. 


Note: ¢ is the quotient (distance from level of widest part of thermal to virtual origin) 
+ (distance from front to virtual orgin). 


It may be shown that the vortex-ring model leads to the equation 
(Mg/p)? t. (8) 

This is consistent with (6) and (7) ifa~3bm = (0-737)? for all thermals. It is natural 
to suppose that the value of m, in each experiment is related to the manner of 
release. The controllable variables are the initial volume of the buoyant mass, 
the radius of the cup, and the speed with which the cup is overturned. The same 
cup (radius 2in.) was used throughout. The density of the solution released was 
roughly 1-3 g/cm? in most experiments, so that the initial volume was roughly 
proportional to the mass excess for the experiment. However, according to the 
set of experimental results presented in table 1, there is no evidence of a relation- 
ship between the configuration at release and the corresponding value of 7. 
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Preliminary trials showed that if the cup was inverted slowly, part of the fluid 
mass moved appreciably while part was still retained by the cup. The release then 
had the character of a pouring, and this tended to initiate plume-like rather than 
thermal motion. The release was so rapid in the main experiments that there was 
no visible distortion or descent of the buoyant mass while the cup was being 
overturned. However, the rate of overturning was not closely controlled or 
measured, and it is therefore possible that variations in this rate may have 
affected the value of n,. In either case we should conclude that a wide variety of 
thermal motions may result from only slightly different initial circumstances. 


3. Separation at the interface 


The junction between the two layers of surrounding fluid was not strictly a 
discontinuity of density. However, we shall assume that the effect of the interface 
upon the thermal is determined only by the over-all difference in density between 
the layers, Ap. This assumption implies that, whenever any part of the thermal 
was being appreciably affected by the interface, a vertical dimension of the 
thermal should be large compared with the depth over which the density of the 
surroundings varied. This thermal thickness was at least four times the interface 
thickness in these experiments, so that the assumption seems reasonable. 

If V, is the volume of the thermal as its front leaves the upper layer, the mean 
density excess in the thermal is then J//),. Material having this density excess 
would be neutrally buoyant in the lower layer if W/V, = Ap, that isifV, Ap/M =1. 
An increase in the parameter V; Ap// therefore increases the difficulty with which 
a thermal may progress into the lower layer. If V, Ap/J/ is sufficiently great, the 
whole thermal will be buoyant in the lower layer, and must eventually return 
towards the interface. The magnitude of this critical value of V, Ap// will depend 
on the rate at which the particular thermal entrains external fluid during its 
passage through the interface, and on the distribution of the intensity of this 
entrainment over the thermal surface. Some allowance for these factors will be 
obtained if the parameter is replaced by VAp/J/, where V is the volume of the 
thermal when its widest part is just leaving the upper layer. Let 


VAp|M = 8S. (9) 


Also let NV, be the number of marked particles released, and let N, be the number 
which would still be moving downwards as t-> 00 if the lower layer were of 
unlimited depth. Let NIN, = Y, (10) 


where Y will be called ‘the yield’. Evidently Y = 1 when S = 0, since the whole 
motion then takes place in unstratified surroundings. By the argument already 
given, Y = 0 for S > £, where £ is a positive constant. We therefore distinguish 
two classes of thermals: 


‘strong’ thermals, for whichO S<f and Y= (11) 


‘weak’ thermals, for which S>f and Y=0. 


The experiments were intended to determine the function F and the value of /. 
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Experimental method 


The substance released was a solution of cobalt chloride in water, so that in our 
case N, was the number of cobalt ions released. The solution was strongly coloured, 
and visual and photographic observations were therefore possible without the use 
of added dyes. 

After each experiment, the cobalt had become divided into two layers. The 
first layer was relatively shallow and was situated about the interface. This layer 
therefore contained (NV, —N,) ions of cobalt. The second, deeper, layer rested on 
the bottom of the experimental tank, about 1-5 m below the interface. This layer, 
which contained the ions which would have penetrated indefinitely into an 
unbounded medium, was discarded by partly draining the tank. The remaining 
contents of the tank were stirred thoroughly so as to produce a homogeneous 
solution, the volume of this solution (about 3201.) was found geometrically, and a 
number of 20 ml. samples were removed. 

The concentration of cobalt in the samples was estimated by the method 
originally proposed by Haywood & Wood (1943) as modified by Eborall (1946). 
Briefly, the cobalt ions in the solution were caused to form a chemical complex 
by reaction with ‘sodium nitroso-R salt’ (sodium 1 nitroso 2 hydroxynaph- 
thalene-3:6 disulphonate). The resulting substance, in weak solution, has an 
absorption coefficient for blue light which is directly proportional to its concen- 
tration. A measurement of the absorption over a known path length therefore 
gave the concentration of cobalt in the samples, and so the total amount of cobalt 
which had collected near the interface was determined. The method is not 
absolute, and the necessary calibration was carried out by mixing known small 
volumes of the stock solution of cobalt chloride into 3201. of water, sampling and 
estimating the cobalt in the manner just described. The amount of cobalt re- 
maining near the interface after an experiment was then expressed as a volume 
of the stock solution to which it was equivalent. This volume was divided by the 
volume of solution released to obtain the value of (1— Y) for the experiment. 


Determination of Ap 


The difference in density between the two layers was usually about 10-° g/em°. 
A long uniform glass tube was partly filled with paraffin wax and weighted so as 
to weigh about 20 g when sealed at both ends and submerged in water. This sinker 
was suspended in the tank by two nylon monofilaments, one of which, attached 
near one end of the tube, was used to make the tube axis horizontal, as indicated 
by a spirit level fixed inside the tube. The other filament was attached near the 
middle of the tube and led to one arm of a balance. An application of the principle 
of moments determined what proportion of any change in the upthrust on the 
sinker would act on the filament attached to the balance. The tank contents being 
at rest, a change in the weight of the sinker was related, by Archimedes’ principle, 
to a corresponding change in the ‘mean density’ of the surroundings. Terms of 
odd order in the distribution of density with depth did not affect the upthrust 
because the axis of the sinker was horizontal and the sinker was geometrically 
symmetrical about its axis. If any terms of even order are sufficiently small, the 
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mean density found from the upthrust approximates to the density at the level 
of the axis. This approximation was always sufficiently accurate in these 
experiments. 

The sensitivity of the arrangement was 1g change in balance for about 
10-%g/em* change in density. The sinker was accordingly weighed at various 
depths, a typical graph of density change against depth being shown in figure 3. 
It will be seen that the interface was about 3 cm thick in this case. 
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Figure 3. An example of the variation of the density of the surrounding fluid with depth. 
In the upper layer (low values of z) the density was constant. 'The lower layer had a slightly 
stable stratification. 


Determination of V 


Two cameras, so placed that their optic axes were respectively perpendicular to 
the front wall and to one of the side walls of the tank, made simultaneous exposures 
at the instant when the widest part of the thermal profile seen from the front of 
the tank was level with the top of the interface. Profile sketches like those shown 
in figure 4a were made from the enlarged images of these photographs. The 
profiles often indicated marked departures from axial symmetry. 

With respect to the axes shown in figure 4a, Av and Ay were measured at 
several levels z, in the thermal. Then V was found from the approximate expression 


[ae Ay) (12) 


The right-hand side of (12) was evaluated as 17 times the area under a graph such 
as figure 4b. Equation (12) is accurate in the case of axial symmetry, and may be 
accurate, depending on the orientation of axes, when horizontal cross-sections of 
the thermal are elliptical. 
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Experimental results 


The graph of Y against S is given in figure 5, which shows that Y varies roughly 
linearly with S over a certain range. The slope of the line is approximately — 0-5. 


The results in figure 5 and table 2 may be expressed in the form 
Y=0-95-438 when 01<S< a3 
when S > 1-9. 


Additionally, we know that Y = 1 at S=0, but the experiments have not 
revealed the relation between Y and S in the range 0 < S < 0-1. 
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FicureE 4. (a) Simultaneous profiles of a typical thermal. (b) Graphical determination of 
the thermal volume V (Aa and Ay are the widths across the profiles at the level z,). 


Now, S is not the only dimensionless parameter which can be constructed from 
the conditions of the experiments. For example, h/V} is another such parameter, 
where h is the vertical distance of the interface from the virtual origin of the 
thermal motion. Although h/V} varied between 1-0 and 4-1 in these experiments, 
the variations appear to have had no significant effect upon Y. 


4. Motion in the lower layer: weak thermals 

In this case there is a lower limit tothe progress of any part of the thermal mass. 
Outline sketches of a typical weak thermal at successive instants are shown in 
figure 6. The distance z of the front from the virtual origin of motion in the upper 
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layer increases to a maximum (at the ‘culminating point’) and then decreases. 
Let s be the value of z, and t, the value of t, at the culminating point. It was found 
that, for every weak thermal when t < ¢,, the graph of (z—s) against (t—t,)? was 
a straight line in the lower layer. An example is shown in figure 7. 

This was motion with a constant deceleration, such as might be expected when 
a parcel of fluid moves through uniform surroundings in the direction opposed to 
that of the buoyancy force. Since mixing at the front between the parcel and its 


104 


y 054 


Ficure 5. Relation between the fraction of the mass released which penetrated the 
inversion and the relative strength of the thermal. 


Ap 
M (g/c.e. V k, h a, 
Expt. (g) x (l.)  (em%/sec) 7, (em) (em) (em/sec?) v/v, 
6 8-07 5-62 3:79 594 4-2 16-5 31 0-66 0-48 
17 11-4 5-89 4-01 505 3:8 31:3 33 0-32 0-59 


51 1-05 1-64 5-72 224 4-4 8:8 39 0-15 0-58 
52 2-10 1-64 2-84 572 6-7 16-4 57 0-21 0-53 
53 3-10 1-62 3-78 223 3-1 25:3 29 0-14 0-73 
54 2-48 1-62 4-98 154 3-1 16-6 30 0-17 0-92 
55 0-64 1-69 4:36 74 2-6 4:7 24 0-14 0-73 
56 1-48 1-69 3°81 304 5-0 17-2 38 0-29 0-80 
57 1-25 4-12 2-20 248 51 7:9 36 0-42 0-74 
58 2-89 4-12 2-30 — 12-9 0-66 

59 4-08 4-01 1-81 624 5-9 27-2 39 0-23 0-60 
60 6-66 4-01 4-32 121 2-2 20-9 16 0-58 1-28 
61 4-16 4-29 4-29 138 21 9-9 20 0-69 1-06 
62 6-22 4-29 3°98 760 5:5 22-9 46 0-35 0-48 


TABLE 2. Measurements on weak thermals (Y = 0). 


Note: x is the distance from the interface to the culminating point. 
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FicurE 6. The motion of a weak thermal above and below an interface. The solid curve 
shows the (z, t)-variation, and the other symbols distinguish the successive chermal profiles. 
Low values of z refer to the upper layer. 
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Figure 7. The constant deceleration of the front of a weak thermal 
while in the lower layer. 
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surroundings then has the effect of working more rapidly against the buoyancy 
force, the intensity of such mixture may be expected to decrease rapidly as the 
thermal moves into the lower layer. 

In the lower layer, oem: (14) 


Accordingly, let 2 be the distance from the interface to the culminating point. 
If (14) held near the interface, the velocity of the front at the interface would be 


Ve = \{(—2a,2), (15) 


which will be called the ‘virtual velocity’. But, from (6), the actual velocity of 
the front at the interface is o, = k,/2h, (16) 


since z = h at the interface. 

It will now be assumed that the particle of fluid which is at the front when the 
weak thermal touches the interface remains at the front until the culminating 
point is reached, and that there is no mixture between other fluid and this particle 
while it is in the lower layer. 

The density difference between the particle and its environment, which consists 
largely of fluid of the lower layer, is 


(K, M/V)—Ap, (17) 


where A, is a numerical constant. 

The kinetic energy per unit volume of the particle is some fraction of the 
average potential energy per unit volume given up by the thermal in its progress 
from the virtual origin to the interface. The kinetic energy per unit volume of the 
particle at the interface is therefore 


0, ViMg/V, (18) 


where C, is a numerical constant. If the energy which the particle gains while in 
the lower layer is negligible, the kinetic energy density (18) is equal to the 
potential energy density of the particle at the culminating point, so that 


_{(K, M|V)—Ap} ga = C,V3Mg/V. 
That is, Vi/a = (S—K,)/C,. (19) 


Experimental results 


These are shown in table 2. The value of a, was determined as the slope of the 
corresponding graph like figure 7, x was measured directly from the photographic 
images, and v, and v, were calculated from (15) and (16). 

It will be noticed that v,/v, varied between 0-48 and 1-28. This seems consistent 
with Miss Woodward’s observation (1959) that the maximum mean vertical 
component of the fluid velocity in a thermal moving through unstratified sur- 
roundings is about twice the velocity of the thermal front. 

The experimental values of V/x have been plotted against S in figure 8, for both 
strong and weak thermals. It will be noticed that there is a sharp break, corre- 
sponding to the division between strong and weak thermals, near S = 1-9. In the 
weak-thermal region, to the right of this break, a roughly linear relationship is 
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evident. If the line were produced, it would pass through or near the origin, and 
V}/ais therefore roughly directly proportional to S. This is the form taken by (19) 
if A, is much less than the minimum value of S which is possible in weak thermals, 
namely S = 1-9. The value of C,, obtained from the slope of the graph in figure 8, 
is 3-6. 


Ficure 8. The penetration of an interface by thermals. Points to the left and right 
of the step represent strong and weak thermals respectively. 


5. Motion in the lower layer: strong thermals 

Substantially all that part of a strong thermal which would proceed indefinitely 
into the lower layer is denser than the surrounding fluid. If the advancing part 
had not been stopped, the momentum which it had on leaving the interface would 
eventually have become a small fraction of the total momentum, which would 
have been continually augmented by the action of the buoyancy force. It there- 
fore seemed possible to suppose that thermal motion might be resumed in the 
lower layer at a sufficient distance below the interface. 

This motion would be described by the equations 


Z = No", (20) 
where k, = (M,e/p,)}. (21) 


Equations (2), (6) and (7) have been used with a suffix 2 to refer to variables in the 
lower layer. C, would have the value obtained from thv experiments in the upper 
layer. 

Experiments 


The observations in each experiment were reduced to graphs of (z —z,) against 2r 
and of z? against ¢. Typical results are shown in figures 9 and 10. In both these 
experiments it will be seen that z = n,7 throughout the motion in the lower layer, 
verifying (20). In most experiments, however, r seemed to be nearly constant in 
the earliest part of the motion in the lower layer, while the advancing part was 
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separating from the remainder of the marked fluid. The first equation of (21) was ny) 
always approximately obeyed in the lower layer, typical cases being shown in b 
figures 9b and 105. It is therefore reasonable to suppose that (20) was also always e: 
approximately obeyed, but that the lateral extent of the advancing part was 
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obscured by the remainder for a varying interval. In some cases this transitional 
behaviour persisted throughout the observed motion in the lower layer, so that no 
estimate of n, could be obtained. 

From (6), (7) and (21), if thermal motion is resumed in the lower layer, we have 


kynj/kyng = 


since p = p, inall experiments. The mass excess in the lower layer was always less 
than J since, first, the yield was always less than unity, and secondly, the 
surrounding fluid in the lower layer was always denser than that in the upper, so 
that the mass displaced by the advancing part of the thermal was increased and 
its mass excess therefore reduced. Hence, if thermal motion is resumed in the 
lower layer, 3 

4 kant |kyn <1. (23) 
The experimental values of k,n?/k,n§ are shown in table 3, where it will be seen 
that about half these values are greater than unity. 


ky ky 
Expt. (em?/sec) No kyni/kyng Expt. (cm?/sec) No kyni/kyn} 

a 584 3-6 0-81 26 283 Ee? 1-41 

8 680 4-1 0-79 27 708 — — 

9 531 28 459 3:8 1-00 
1l 958 4:7 1-15 29 319 — — 
12 1238 3-4 1-23 30 382 3-1 1-17 
13 704 — — 31 447 4-2 1-19 
14 710 4-2 0-99 32 613 34 1-15 
16 551 3:0 1-44 34 576 4-6 0-65 
21 720 35 282 2-3 1-13 
22 499 3-1 0-92 36 412 3-4 0-93 
23 834 4-1 1-15 38 388 3-2 0-93 
24 527 2-7 1-31 39 270 2-1 1-68 
25 331 2-3 0-87 


TABLE 3. Comparison of motions above and below an interface (strong thermals). 


We conclude that the motion of the advancing part of a strong thermal differs 
from thermal motion from rest at a horizontal surface into unstratified fluid. 
Returning to the vortex-ring model, we see from (8) that values of k,n} /kyn3 
greater than unity may occur if, and only if, the value of a~*bm is sufficiently 
increased by the passage through the interface. This means that, according to the 
vortex-ring model, the distributions of velocity are sometimes necessarily 
different in the two layers. 
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for this research and for permission to publish the results. 


: 

— 


384 


J. M. Richards 


REFERENCES 


BJERKNES, J. 1938 Quart. J. Roy. Met. Soc. 64, 325. 

Davies, R. M. & Taytor, G. I. 1950 Proc. Roy. Soc. A, 200, 375. 
Esporatt, F. E. 1946 Metallurgia, 39, 104. 

Haywoop, F. W. & Woop, A. A. R. 1943 J. Soc. Chem. Ind. 62, 37. 
Lams, H. 1932 Hydrodynamics, 6th ed. Cambridge University Press. 
Scorer, R. 8. 1957 J. Fluid Mech. 2, 583. 

Scorer, R. 8. 1958 Natural Aerodynamics. London: Pergamon. 


Scorer, R. 8. & Lupiam, F. H. 1953 Quart. J. Roy. Met. Soc. 79, 94. 


StomMEL, H. 1947 J. Met. 4, 91. 
TURNER, J. S. 1960 J. Fluid Mech. 7, 419. 
Woopwarp, B. 1959 Quart. J. Roy. Met. Soc. 85, 144. 


|_| 

ec 
de 
flc 
to 
pe 
ve 
Ww 
he 
flc 
by 
th 
th 
| W 
| de 
at 
co 
an 
4 dis 
de 
sic 
(1 
a wi 
th 
fla 
pr 
for 
th 
pa 
fro 

| 


385 


The effect of a sudden density change on a slightly 
non-uniform flow 


By W. R. HAWTHORNE Anp P. J. BANKS 


Engineering Department, University of Cambridge 
(Received 22 June 1960 and in revised form 24 March 1961) 


The effect of combustion or heating on a slightly non-uniform flow in a tube of 
constant radius is explored by examining a simplified model, in which a density 
decrease occurs across a plane actuator-disk normal to the axis of the tube. The 
flow is considered to be steady, axi-symmetric, incompressible and inviscid and 
to originate at an injector-plate which is situated upstream of the disk. This 
paper analyses and discusses the effects caused by the density decrease when the 
velocity of flow through the injector-plate has a fractional non-uniformity w,, 
when there is a curvature @, of the injector-plate, and when the density decrease 
has a fractional non-uniformity 6. It is shown that the non-uniformity of the 
flow far downstream of the disk caused by the disturbances wu, and 9, is diminished 
by the density decrease; while that due to the disturbance 0, is increased, unless 
the density decrease occurs close to the injector-plate. The distance separating 
the injector-plate and actuator-disk is found to be an important parameter. 
When this distance is small, compared to the radius of the tube, and the density 
decrease is severe, a small disturbance w,, 0, or 6, may set up a pressure variation 
at the injector-plate of many times the inlet dynamic pressure of the flow. It is 
conjectured that this pressure variation could affect the uniformity of velocity 
and composition of the flow through the injector-plate, and thus cause the 
disturbances u, and 6, to grow or diminish. 


1. Introduction 


In rockets, ramjets and turbo-jet combustion chambers a flow with large 
density changes occurs. The decrease in density caused by combustion con- 
siderably affects the pattern of flow. Such effects have been studied by Scurlock 
(1948), Tsien (1951), Emmons, Ball & Maier (1952) and others, in connexion 
with flow through a V-shaped flame-front stabilized at the centre of a duct. In 
these studies, the effect of the flame-front on the flow is due to the shape of the 
flame-front rather than to any non-uniformity of the approaching flow.+ In 
practical cases the flow approaching the flame-front is generally non-uniform; 
for instance in rocket combustion chambers the flow of oxidant and fuel 
through the injector-plate varies with radial distance. 

The effect of combustion or heating on non-uniform flows is explored in this 
paper by examining a simplified model in which a density change occurs across 

t Ball (1951) considered the effect of shear in the flow approaching a V-shaped flame- 
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a plane actuator-disk normal to the main-flow direction. This actuator-disk 
would represent a flame-front if the non-uniformity of the flow at the front was 
identical with the flame-speed variation in the approaching flow, caused for 
instance by non-uniform mixture strength. 

In the model examined (see figure 1), a steady, axi-symmetric, incompressible, 
inviscid and slightly non-uniform flow in a tube of constant radius approaches 
a density-change actuator-disk normal to the axis of the tube. On the upstream 
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FicgureE 1. Diagram showing flow model and nomenclature. 


side of the disk the flow is assumed to be of uniform density. Downstream the 
density may be either uniform or may vary owing to a slight non-uniformity in 
the density change at the disk. In accordance with the assumption of inviscid 
flow there will be no heat transfer, so that, since the flow is also assumed to be 
incompressible, the density may be considered to be constant along a streamline, 
except across the disk. Thus the flow on either side of the disk is obtained by a 
small perturbation of a flow with uniform velocity and density, and the problew 
may be linearized by neglecting the products of perturbations in the equations 
defining the flow. The effect of gravity on the flow is neglected. The disk is sup- 
posed to be stabilized at a finite distance downstream of a plate through which 
the fluid is injected, and the flow is assumed to continue for an infinite distance 
downstream of the disk. 


2. Solution for flow on either side of the density-change actuator-disk 


The type of flow, downstream of the disk, described above is termed stratified 
or non-homogeneous; the flow upstream may be regarded as a particular case of 
such a flow. Long (1953) has shown that the vorticity 7 in a steady, stratified, 
plane two-dimensional flow varies along the streamlines as the velocity, of 
magnitude q, changes, in accord with a relation of the form 
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where p is the density, ‘ a stream function such that d'¥/dn = —q with n 
measured normal to a streamline, and H is constant along a streamline. Equation 
(1) is equation (8) of Long (1953) with 7 replacing (— V?P), q? replacing (V‘Y’)?, 
and the term due to gravity omitted. The equation is also applicable to axi- 
symmetric flow with the stream function defined by d'Y/dn = —qr, because of 
the similarity between axi-symmetric and plane two-dimensional flow, The 
vorticity 7 is then the circumferential component. 

In the problem considered in this paper, the flow is axi-symmetric and involves 
a small perturbation of a uniform flow, with velocity U and density p, in the 
axial direction. We define the velocity components in the axial, x, and radial, r, 
directions to be U(1+u) and Uv respectively, and the density to be p(1+94), 
where uw, v and 6 <1; then, by substituting in equation (1) and neglecting 
products of perturbations, one obtains the equation 


dé 
y= H+ av (2) 
Now JU is constant and H and 6 depend only on ‘’, so that the vorticity is con- 
stant along the streamlines of the perturbed flow, to first order in the perturba- 
tions. These streamlines depart only slightly from the streamlines of the un- 
perturbed flow, which are lines of constant radius. Hence in the perturbed flow 
the vorticity is a function of radius alone to first order in the perturbations, and 


may be expressed as 
Cv Ou, 
y= (3) 


where Uu,, is the value of the axial velocity perturbation at a plane where v and 
dv/ox = 0. 

The flow may be considered to be the sum of the unperturbed flow and a per- 
turbation flow, and, since the products of perturbations are being neglected, the 
perturbation flow may be resolved into a number of parts. Thus if we write 
u = U.+u', where w’ and v both approach zero at the same plane, then equation 
(3) shows that Uw’ and Uv are the velocity components of the irrotational part 
of the perturbation flow. The rotational part is seen to have axial velocity Uw,, 
and zero radial velocity. Therefore a perturbation stream function y, defined by 


ru’ and rv = 


satisfies Laplace’s equation, which in axi-symmetric co-ordinates is 


er\r or) dx? 


The solution of this equation which satisfies the boundary conditions, v = 0 at 
r = 0 and r = 1, the radius of the tube containing the flow being taken as the 
unit of length, is 
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where A, and B, are arbitrary constants, and k,, are a series of constants deter- 
mined from 
= 0, (7) 


so that k, = 3-832, k, = 7-016, k, = 10-174, ete. Therefore the perturbation flow 
may be considered as the sum of a series of Fourier—Bessel harmonics. The nth 
harmonic of the perturbation stream-function, and of the radial velocity, has 
a radial form defined by the Bessel function J,(k,,7), which has (n—1) zeros 
between the tube axis and the tube wall. The factor (in any harmonic) which is 
independent of the radius r, e.g. ky 1{A, B, in equation (6), will be 
referred to as the amplitude of the harmonic. 

Using equations (4) and (6) the velocity perturbations Uw and Uv in the flow 
on either side of the density-change actuator-disk are found to be given by expres- 
sions of the form 


w= Unk, r); (8) 
n=1 
and v= > (9) 
n=1 
where the amplitudes wu, and v,, of the nth harmonics are 
Un With uy = —(A, +B, e-kn®) (10) 
and U, = (A, — B, (11) 


The amplitude w,,,, of the nth harmonic of u,, is defined by equation (8). The 
quantities wv and v will be referred to as the axial and radial velocity perturbation 
fractions respectively, and w will be used to define the non-uniformity of axial 
velocity. 

An expression for the pressure-perturbation fraction p, that is the static 
pressure-perturbation divided by the dynamic pressure of the unperturbed flow, 
may be derived by equating the terms containing perturbations on either side 
of the equation of motion in the radial direction, and neglecting products of 
perturbations, from which one obtains 


—(4pU2) = (pU 


(12) 


from which it may be concluded that 
(13) 


the constant vanishing because, at the plane where w’ = 0, the flow is parallel 
to the axis. 


3. Application of boundary conditions 


The quantities w,,,,, A, and B, in the expressions for the perturbation fractions 
are determined by the boundary conditions at the injector-plate, at the density- 
change actuator-disk, and at an infinite distance downstream. The additional 
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suffixes 1, 2, a, and e will be used to denote conditions upstream of the actuator- 
disk, downstream of the actuator-disk, at the actuator-disk, and at the injector- 
plate respectively. 

The boundary conditions at the injector-plate are determined by the non- 
uniformity of the flow through the plate and the geometry of the plate. The plate 
is assumed to be symmetrical about the axis of the tube, and to be slightly curved. 
Hence the perturbation of the axial velocity at the plate U,u, arises directly 
from the non-uniformity of the flow through the plate, and can be defined by 


UenI (14) 

n=1 
The curvature of the plate will produce a radial velocity at the plate, which, as 
a fraction of U,, may be written 


= 0, = A key”), (15) 
1 


where @, is the angle by which the entering flow is deflected inwards from the 
axial direction due to the curvature of the plate. It is equal to the angle by which 
the surface of the plate deviates from that of a flat plate. Near the tube wall 
§, must approach zero so that the flow adjacent to the wall is injected parallel 
to the wall. 

The flows upstream and downstream of the density-change actuator-disk are 
matched at the disk by the conditions of conservation of mass and momentum. 
From conservation of mass, we derive that 


U,(1 + = Po(1 + 63) + Ugg). 


By neglecting the term containing a product of perturbations, and equating the 
terms corresponding to the unperturbed and the perturbation flows on either 
side of the equation, one obtains the results 


PLU, = (16) 
and Uay = Ug2t (17) 


Hence, using equation (8), the non-uniformity of the density change may be 
expanded into harmonics, to give 


= Bon 7). (18) 


From the conservation, across the disk, of mass and radial momentum, it is 
apparent that the radial velocity is continuous across the disk, so that, using 


equation (16), 
Nar = Var (19) 


where A = p,/p, defines the ratio of the densities of the unperturbed flow on either 
side of the actuator-disk. 

Conservation of momentum normal to the plane of the disk gives an equation 
which, on neglecting terms containing products of perturbations, equating the 
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terms corresponding to the perturbation flow on either side of the equation, and 
using equation (16), becomes 


MDa 2(2Uge + = 2AUq1); 
which, with equations (13) and (17), reduces to 
— Ung = (1L—A) Ug. (20) 


This equation defines the change, at the actuator-disk, of the rotational part of 
the perturbation flow, w,,, which is otherwise invariant along the tube, in terms 
of the conditions at the disk. Thus it leads directly to an expression for the change 
in the vorticity of the flow at the disk, and the relation of this result to the work 
of other authors is discussed in § 4.2. 

At an infinite distance downstream from the disk, the streamlines are assumed 
to be parallel. Thus at 2 = 00, u’, vand 0v/dx = 0, and U,u,,2 is the only remaining 
velocity perturbation, so that, from equation (11), 


(21) 


By substitution of the solutions for the velocity perturbation fractions, equa- 
tions (8) to (11), in the boundary conditions, equations (14), (15), (17), (19) and 
(20), and by use of equations (18) and (21), expressions for w,,,, A, and B,, can 
be obtained in terms of A, e~*n4, u,,, O,, and 6,,, where a is the distance separating 
the injector-plate and actuator-disk. These expressions are shown in the Appen- 
dix, in which w,,,, A, and B,, are related to u,,, 0,,, and 6, by influence coeffi- 
cients. The equations stating the dependence of wu, v and p on 2, r, Un», A, and B, 
are also summarized in the Appendix. Using the Appendix, the flow pattern 
arising from any combination of the disturbing non-uniformities, w,, 0, and 6,, 
can be assessed. The number of terms required in the summations for wu, v and p 
will depend on the number of Fourier-Bessel harmonics required to form each 
disturbing non-uniformity. 


4. Discussion of solution 
4.1. Assumption of incompressible flow 


The effect of the compressibility of the fluid will become appreciable when the 
Mach number of the unperturbed flow downstream of the density-change 
actuator-disk exceeds a value of about 0-3. The ratio of the Mach numbers of 
the unperturbed flow on either side of the disk is given by 


M 
if the fluid is assumed to be a semi-perfect gas, since then p, 7, = p27, for small 


M, and M,. Hence the solution is adequate for 


A < (22) 


If the temperature downstream of the density change were 4500°K, a very 
high value for flame temperature, then A would equal 0-07 and the above con- 
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dition would require M, < 0-08. Combustion reactions in gas mixtures proceed 
at Mach numbers much less than this value, so that the assumption of incom- 
pressible flow is justified for the flow of gas mixtures normally through an 
actuator-disk at which combustion takes place. 

In rocket combustion chambers, the Mach number of the flow when combustion 
is completed is determined by the ratio of the cross-sectional areas of the com- 
bustion chamber and the nozzle throat, where M = 1. If this area ratio equals 
3 or more then M, < 0-3 and compressibility effects may be neglected. 


4.2. Vorticity change at density-change actuator-disk 


The change in the vorticity of the flow at the disk is given by equation (20) with 
equations (3) and (16), from which 


(1 A) 


(23) 


Hence the vorticity of the flow is changed at the disk by an amount determined 
by the radial gradient of axial velocity just upstream of the disk and the density 
ratio across the disk. The vorticity change is zero if the density change is small, 
and does not depend directly on any small non-uniformity in the density change 
at the disk, to first order in the perturbations. 

Hayes (1957) has studied the general case of rotational compressible flow 
through a surface of discontinuity in density. Hayes shows that only the com- 
ponent of vorticity tangential to the surface experiences a change as the flow 
passes through the surface, and that this change is made up of two parts, both 
of which depend on the ratio of the densities on either side of the surface. One 
part of the change is proportional to the gradient along the surface of the mass 
velocity normal to the surface, and this part reduces to equation (23) for the 
conditions of the problem considered in this paper. The other part varies with the 
product of the component of velocity tangential to the surface and its gradient 
in its own direction, which in this paper is v(0v/ér) and so is zero to first order. 
This latter part provided a vorticity change at the flame-fronts studied by the 
authors cited in $1. In the papers by these authors, the flame speed, which is 
the velocity of flow normal to and into a stationary flame-front, and the density 
upstream of the flame, were considered constant, so that the part of the vorticity 
change which is important in this paper was identically zero. 


4.3. Case of injector-plate and actuator-disk at an infinite 
distance apart (a = 0) 
If the injector-plate is at an infinite distance from the actuator-disk, then there 
is some point between them at which the radial velocity approaches zero. Hence 
the relevant disturbing non-uniformities in this case are upstream vorticity, 
represented by w,,;, and non-uniform density change at the disk, d,. The upstream 
vorticity may be caused by either or both of the disturbances u, and @, at the 
injector-plate, and the expression for the amplitude of the nth harmonic of w,., 
is, from the Appendix, 


Ucin = Ben: (24) 
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The main results of the analysis in this case do not depend on the form of radial 
variation of the disturbances. For with a = 0, so that ¢€, = e~*»¢ = 0, it is seen 
from the table in the Appendix that w,,,, A, and B, are independent of k,,. 
Therefore at x = —0oo,x = Oand x = +, where ekn? = 0, etn = Land e-*knt = 0 
respectively, the amplitude of the velocity perturbation caused by a disturbance 
consisting of a single harmonic of unit amplitude is the same for all harmonics. 
Thus the radial form of a velocity perturbation is identical at the disk and at an 
infinite distance upstream and downstream of the disk, even if the disturbance 
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Ficure 2. Axial velocity perturbations caused when the injector-plate is at an infinite 
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consists of more than one harmonic. For example the amplitude of the nth 
harmonic of the axial velocity perturbation fraction just upstream of the 
actuator-disk is, from equation (10), 

Uain = ~ Ary ~ Big, 
which, using the influence coefficients in the Appendix, becomes, with e, = 0, 


Uan = (ten Cun) (1 A) (ten + Ben) + 


and, with equations (24), (8) and (18), this can be summed to give 


Ugqy = + 463. (25) 
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The variation of the axial velocity perturbation fraction, expressed as u/u,.,, 
and of the perturbation, expressed as Uu/U,u,,;, in the axial, x, direction along 
a line of constant radius takes the form shown in figure 2, which shows separately 
the effects of w,,, figure 2(a), and 6,, figure 2(b). The perturbation fraction is 
indicated by full lines, and the perturbation by full lines upstream of the disk 
and broken lines downstream. The ordinates of significant points on the curves 
are shown in terms of the density ratio A, and the curves are drawn to scale for 
A = 0-4. The scale of x is omitted since the various Fourier—Bessel harmonics 
of a perturbation differ in their manner of variation with x, though the harmonics 
behave identically at 2 = —o, x = 0 and # = + as explained above. 

At equal distances on either side of the disk the irrotational part of the per- 
turbation Uw’ is equal in magnitude and opposite in sign, for either of the dis- 
turbances 6, and w,,,. From figure 2 it may be seen that these disturbances cause 
identical values of Uw’ at the disk, and thus everywhere in the flow, for a par- 
ticular Fourier—Bessel harmonic, if 


6, = —2(1-A) wa}. (26) 


Identical values of Uw’ lead to identical values of v and p, from equation (12). 
Hence the radial velocities and pressure perturbations in the flow produced by a 
disturbance u,,, may be identically reproduced by a disturbance 6, which satisfies 
equation (26), and vice versa. The exceptional case A = 1 is discussed in § 4.4.4 
in connexion with the similar result obtained when a is finite. 

From figure 2(a), for which u,, + 0 and 6, = 0, it is seen that a disturbing 
non-uniformity of axial velocity w,., is severely attenuated by density decrease, 
SINCE Uy» = A®u,,; the axial velocity perturbation is also attenuated by the 
ratio A. The axial variation of the perturbation is seen to be of very simple form, 
the perturbation just upstream of the disk being equal to that far downstream, 
and vice versa. 

From figure 2 (b), for which u,, = 0 and 6, + 0, it is seen that a disturbance 64, 
representing a non-uniform density change at the disk, causes an irrotational 
upstream flow, w,., = 0, to become rotational downstream of the disk. The magni- 
tude of the rotational part of the downstream axial velocity perturbation 
fraction or non-uniformity w,,. increases as the density ratio A decreases, and 
approaches the maximum value of (36,) as the downstream density approaches 
zero, A > 0. 


4.4, General case 
4.4.1. The effect on the results of the radial form of a disturbance 


The distance a separating the injector-plate and the density-change actu- 
ator-disk is an important parameter in the problem, and appears in the results 
of the analysis in the quantity ¢,, = e~*x“. Hence the effect of a on the velocity 
perturbations caused by a disturbance consisting of a particular Fourier—Bessel 
harmonic depends on the number n of the harmonic. If a disturbance con- 
sists of more than one harmonic, then, in general, the radial form of a resulting 
velocity perturbation will differ at the injector-plate, the disk, and infinity down- 
stream, in contrast to the result obtained in the special case a = 00, § 4.3. The 


4 
| 
| 
ial 
en 
n° te 
= 
0 
ice 
CS. 
an 
3 
| 
a 
cee 
= 
= 
4, 
4 
: 
we 
| 


394 W. R. Hawthorne and P. J. Banks 


factor k,,a, which determines the effect of a on the results for the nth harmonic 
of a disturbance, is approximately equal to 7 times the ratio of a and the radial 
distance between adjacent zeros of the radial velocity perturbation resulting 
from that harmonic. This rough equality becomes more exact with increase of 
the value of (n—1), the number of zeros of the radial velocity perturbation 
between the tube axis and the tube wall. 

The effect on the results of a being finite rather than infinite is greatest for the 
primary disturbance harmonic, n = 1, so that in the following discussion the 
disturbances are considered to consist of this harmonic alone. The results for a 
disturbance consisting of the nth harmonic alone can be obtained from those for 
the primary harmonic simply by multiplying the values of a by the ratio (k,/k,,). 
The flows resulting from each disturbance acting singly are discussed in the 
following sections. 


4.4.2. Flow caused by non-uniform axial velocity at the injector-plate (disturbance u,) 


The variation of (u/u,) and (Uu/U,u,) with x is shown in figure 3 for values of a 
of 0, $f, B and oo, where # = kj 1In2. The arrangement of this figure is similar 
to that of figure 2(a) and the curves from figure 2 (a), for a = 0, are repeated 
on it. The variation of the axial velocity non-uniformity (w/1,) is shown by the 
full lines, and it is seen that as a is reduced so is the attenuation {1 —(u..2/u,)} of 
the disturbance by the density decrease. The attenuation increases as the 
density ratio A is decreased, for all values of a except zero. In the limiting case of 
a = 0, the density change occurs immediately after injection and the flow is 
invariant along the tube. The axial velocity perturbation (Uu/U,u,) is shown on 
figure 3 by full lines upstream of the disk and broken lines downstream. It is 
seen that the axial velocity perturbation is greater at infinity downstream than 
at the injector-plate for a < #, and less fora > f. 

The variation with a of the upstream value of the rotational part of the axial 
velocity perturbation fraction, w,,,, is also indicated on figure 3. It is seen that 
Ux, Iereases as a decreases; therefore, since wu is nowhere greater than w,, 
U, = (Uy — U1) also increases in magnitude as a decreases. Hence, from equation 
(13), it is seen that the pressure perturbations upstream of the disk increase in 
magnitude as a decreases. The radial pressure variation set up at the injector 
plate by the disturbance is of particular interest, since this variation may alter 
the uniformity of injection, and thus the disturbance, if the rate of injection of 
fluid at the plate is pressure controlled, as in liquid rockets. The pressure per- 
turbation fraction at the injector-plate, p,, is shown on figure 4, plotted against a, 
for various values of A. It is seen that for a < 0-1 and A < 0-2, conditions which 
may occur in liquid rockets, then p, > u,. Hence under these conditions a radial 
pressure variation of many times the inlet dynamic pressure (}p, U7) may be 
obtained at the injector-plate, even if w, is small. 

The maximum value of the radial velocity occurs at the disk, since wu’ defines 
dv/cx, equation (12), and is of constant sign on either side of the disk, changing 
sign at the disk (see figure 3). The variation of v,,, the radial velocity at the disk 
as a fraction of U,, with a and A is shown in figure 5. Because of the difference in 
the radial forms of v,, and u,, the maximum values @,, and “@, are used for com- 
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Figure 5. Radial velocity perturbation fraction at the density-change actuator-disk 
caused by a disturbance u, of primary form. Maximum values are compared. 


4.4.3. Flow caused by a curved injector-plate (disturbance 0,) 


In the absence of density change a curved injector-plate is seen, from equation 
(24), to cause a uniform flow to develop a non-uniformity of axial velocity with 
the same amplitude as 0,. This non-uniformity will be attenuated in the manner 
described in § 4.3 by a density decrease occurring at a plane far downstream of 
the injector-plate. The change in this attenuation when the density decrease takes 
place at a finite distance downstream of the injector-plate is shown by figure 6, 
in which the variation with a of w,9, the non-uniformity of axial velocity far 
downstream of the disk, is shown for various values of A. The attenuation 
{1—(&,,/0,)} of the disturbance is seen to decrease with a except near a = 0, 
and to increase as the density ratio A is reduced, for all values of a. 
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FicurE 6. Non-uniformity of axial velocity far downstream of the density-change 
actuator-disk, caused by a disturbance 9, of primary form. Maximum values are compared. 
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It may be shown that the velocity perturbation fractions ~ and v are every- 
where of the same order of magnitude as the disturbance 6, for all values of a 
and A. However, an appreciable radial pressure variation may be set up at the 
injector-plate, as shown by figure 7 where p, is plotted against a for various 
values of A. It is seen that a radial pressure variation of many times the inlet 
dynamic pressure may be obtained at the injector-plate for A < Q-1 and 
0 <a < 0-1, even if 0, is small. 
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Figure 7. Pressure-perturbation fraction at the injector-plate caused by a 
disturbance @, of primary form. Maximum values are compared. 


4.4.4. Flow caused by a non-uniform density change (disturbance 63) 


It is apparent from the table of influence coefficients in the Appendix that 
the values of A, and B,, and therefore of the irrotational part of a velocity per- 
turbation, caused by a disturbance 6, are identical with those caused by a dis- 


turbance w,, if (27) 


Hence, in view of equation (12), the discussion in § 4.4.2 of the radial velocities 
and pressure perturbations caused by a disturbance u, applies also to the flow 
caused by a disturbance 6,, except when A = 1. The exception arises because if 
there is no density change in the unperturbed flow, the disturbance wu, is invariant 
along the tube, whereas the disturbance 6, necessarily sets up a discontinuity 
in u at the disk (see equation (17)). However, the radial velocities and pressure 
perturbations due to 6, are small when A = 1, so that the conditions for an 
appreciable pressure variation to be set up at the injector-plate by 6, are as stated 
for u,. For density change caused by combustion it is possible that such a pres- 
sure variation may alter the non-uniformity of the mixture ratio of the injected 
fluid and thus cause the non-uniformity of the density change to grow or diminish. 

The non-uniformity of axial velocity far downstream, w,.9, caused by this 
disturbance is shown on figure 8 plotted against a for various values of A. It is 
seen that if there is no density change in the unperturbed flow, A = 1, then u_.5 
varies from zero when a = 00 to the value of — 6, whena = 0. The effect of density 
decrease is to amplify this non-uniformity fora > attenuate it for 0 < a < 
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where # = ky1In2 = 0-181. Thus wu... approaches the value of — 40, forall values 
of a as A approaches zero. 


0 0-1 0:2 


«| 


FicuRE 8. Non-uniformity of axial velocity far downstream of the density-change 
actuator-disk, caused by a disturbance 6, of primary form. 


The authors wish to acknowledge the assistance of Mrs P. Camm in carrying 
out the calculations for the figures. 
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Appendix: summary of solution 
Disturbing non-uniformities 


Ue = Men 0. = Den 
n= n= 


and 
n=1 


Resulting perturbation fractions 
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Influence coefficients 


Uen|E 6.,/E 6,,/E 
r A 1 (1—2A) 
1 1-€, ~ 21—A) 
2 
Bin én 2(1—A) 
Asn 0 0 
rey | A(1 = €n) 
Ban A(1—eén) En té —A) 


Here k, are the roots of J,(k,,) = 0, thus k, = 3-832, k, = 7-016, k, = 10-174, etce., 
A= = and E= {1/(1 2€,, + 
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On corner eddies in plane inviscid shear flow 


By L. E. FRAENKEL 


Aeronautics Department, Imperial College, London 
(Received 11 March 1961) 


Local solutions are found for the inviscid shear flow past an acute corner (of 
included angle < 47 on the side of the fluid) on an otherwise arbitrary boundary. 
Unlike irrotational ‘corner flows’, these solutions are determinate locally, 
provided that the vorticity is known. Under certain circumstances the existence 
of a corner eddy may be inferred. 


1. Introduction 


Yih (1959) has shown that a certain inviscid shear flow in a channel (figure 1a) 
includes closed streamlines. The same phenomenon is seen in the shear flow past 
a semicircular projection on a plane wall (figure 15: details of this solution will be 
given in §4).7 The purpose of this paper is to point out that such eddies may 
always be expected under certain conditions, which will be described, and to give 
simple expressions for the stream function in the immediate neighbourhood of 
the corner. 

We are concerned throughout with plane, inviscid, incompressible shear flow, 
and mainly with the case of constant vorticity. Then, writing 


y= 


for the stream function and coordinates, we have 
= const. =w > 0, say. (1) 


We suppose that at some point the tangents to a solid boundary (or to some other 
streamline) make an angle / < 37, measured on the side of the fluid. Near the 
corner the fluid is sufficiently stagnant for its motion to be determined by the 
rotation, so that (for w > 0) an observer moving with the fluid near the corner 
keeps the boundary on his left. We shall see that the dominant term in the 
expression for the velocity near the corner is independent of the boundary shape 
away from the corner, and that this term dominates no matter what irrotational 
flow past the boundary is added. If now the main flow further out from the 
corner is such that an observer moving with it keeps the boundary on his right, 
a corner eddy occurs. In this case the local solution near the corner implies the 
existence of an eddy, but of course the size of the eddy depends on the flow as a 
whole. 


+ Despite the apparent similarity of figure 1b to figure 1 of Yih (1960) the stream 
functions in question are quite different. 
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Flows with constant vorticity are easy to treat analytically, and viscous flows 
bounded by closed streamlines are known to have constant vorticity in the limit of 
infinite Reynolds number (Batchelor 1956). Accordingly the behaviour of 
solutions of (1) near acute corners will be studied in some detail. 


u 


(= 059) 


(b) 


Figure 1. (a) Yih’s channel flow, (b) shear flow past a semicircular projection. 


We also expect our results to hold for the more general inviscid shear flow 
problem = f(y), (2) 
for if a corner occurs on the streamline i = yo, such that f(y,) + 0, it is pre- 
sumably legitimate to write f(y) ~ f(~%) for the purpose of examining the flow 
very near the corner (where Viv = 0). In the Appendix it is shown that Yih’s 
solution (of an equation of the type (2)) does in fact have the local behaviour 
predicted by this procedure, but no general proof of the occurrence of our corner 
flows in solutions of (2) will be given. 

Dean (1944) has found flows containing eddies in the lee of a peak in certain 
solutions of the Stokes slow-motion equations; such eddies are a different species 
from those considered here. 


2. Corner flows with constant vorticity 

We seek a solution of (1) satisfying yy = const. on 0 = +42 (f < 4m). A parti- 
cular integral is given by y = }wr?, and the harmonic, complementary function 
required by the boundary condition may be found by inspection; thus 


) (2 < 37). (3) 
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The streamlines are the hyperbolae 
cos?(48) sin?(4f) wsin? 
and are shown in figure 2a: an observer moving with the fluid in the acute angles 


keeps the boundary on his left, as was stated above. The solution may be inter- 
preted as a corner flow or as a stagnation-point flow; that is, either AOB or AOC 


Ficure 2. Rotational corner flows. (a) £ < 47; (b) f = 4a. 


could be a solid boundary, or O could be a stagnation point in the body of the fluid. 
The solution (3) is not unique for the region |6| < $/, since conditions have not 
been specified on a closed boundary, and an irrotational flow 


= Cr”? cos (4) 
where C is arbitrary, can be added without violating the boundary condition. 
However, we observe that, since 7/f > 2, 7, dominates for sufficiently small values 
of r. 

The solution (3) becomes infinite for 6 = 47, and we therefore form a linear 
combination of yy, and 7, which remains finite as f > 47. Writing 


po4n cos 


we obtain 
= on +7? logr cos 20-770 sin20) (f=47). 


Again we may add an irrotational flow, yf, with 7// = 2; with C = —o loga/z, 
where a is an arbitrary length, we have 


= wr +7? log (r/a) cos 20 — 70 sin 26} (f=4n). (5) 


The 7? log r term dominates for sufficiently small values of r. 
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The streamlines are shown in figure 2b: they are not drawn for || > 47 because 
7—10y,/00 is discontinuous across the negative real axis (if we choose |(| < 7), so 
that the flow is not a real one there. 

Now let yf be a solution of (1) for any field with a streamline which makes an 
acute angle at some point O, and take coordinates with origin at O and with the 
positive x-axis bisecting the angle. The foregoing considerations make it plausible 
that yy, or the logarithmic part of y, provide the dominant term (apart from a 
constant) in the expansion of yy about 0, in powers and logarithms of r. In the 
next section we shall prove this assertion. 


3. The flow near acute corners on arbitrary boundaries 


Consider a fluid domain Y whose boundary # is a Jordan curve with a corner 
at O of included angle / < 37. Various types of domain are of interest (figure 3): 
here we consider only the case (figure 3a) in which @ is semi-infinite and simply 


¥) 


Ficure 3. Two types of fluid domain to which the present results apply locally. 


connected, since for other cases the proof only requires changes in points of detail. 
The domain can be mapped conformally on to the interior ofa circle (Carathéodory 
1932; Courant 1950) or, alternatively, on to an upper half-plane. Let 


z= = 2(§ +7) 
be one of the latter transformations, such that the corner z = 0 maps on to € = 0, 
and Z@ on to 7 > 0. Let Z be so restricted that we may write, after a suitable 
choice of constants, 
for <3, z= ef? 
y(&, 0) = (E<0, (6) 
= {°sina+O(g) (§>0). 
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A sufficient (but by no means necessary) condition for this is that the two ares of 
B meeting at O have bounded curvature in —d < § < 0Oand 0 < & < 6, respec- 
tively; for then we may first apply the Schwarz—Christoffel mapping appropriate 
to the tangents at the corner, and then the theorem in § 152 of Carathéodory 
(1932). We further restrict Z to have y finite and dy/dx ~ O(«-*) for |x| > « 
then for |¢| 0, 
(k,, k, real). 
Our problem is to solve (1) with the boundary conditions 


y = const. on Z 7 
Uy ow (7) 
~—-wy-U+0, as |z| >oinG. 
v Cy 


C 
Ca 
The solution is 
soy? + (8) 
where y,; is a harmonic function required by the boundary condition on 2%. 
Representing it by sources on 7 = 0, we have 


Fw (G) = log —§); 


w { [84 0 


Here A and B are points on 4 such that cy/cg is continuous for € < €, and 
£ > &,[and O(£-*) for |£| + co], and we confine ¢, to the domain || < ¢, 
0 < arg ¢, < 7, whered > € > 0. For |£| > 6, the integrands in (9) are continuous 
functions of £, uniformly convergent for |g] + 00, and regular functions of &. 
Hence, if F(¢,) denotes a function which is regular in Y,, 
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[2 sin? (a+  A*sin® a + O(A**) 
Now, by the contour integration associated with the beta function, 
= Clog€—(€+ d) log (€+d)+0 (a=1). 
Hence 
Omitting the trivial constant 4 F(0), we have 
2b 
WE, 9) = — sin? (a + +sin? + (0<b<4) 
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Transforming to z, writing z = rexp {i(r + «+ 4/)}, adding the first two terms of 
(8), and simplifying, we obtain 


= lor (1- +000) (0<b<}), (10) 
= on 'r logrcos27+O(r?) (b=4), - (11) 


which is the required result. The sign of yy yields the local flow directions described 
in §§ 1 and 2. 


4, An example 
Consider the flow past a semicircular projection, r = 1, 0 < 0 < 7, ona plane 
wall, 0 = 0, 7, r > 1. The Joukowski transformation 
C=2+21 
maps the fluid domain on to 7 > 0 and the semicircle on to 7 = 0, |&| < 2. For 
the differential equation (1) “nd the boundary conditions (7), we have 


= soy? + U(r—r) sind + Wy, (12) 


wo [2 1 
where p= 0) and = 1-18 
{1 I\?, 
(: log +247]. (13) 


The case shown in figure 1b is for U = 0; adding an irrotational flow with U > 0 
reduces the size of the eddies, and introduces stagnation points on 


1<r<o. 


I am indebted to Miss M. Lawson for computing the streamlines of figures 1b 
and 20. 


Appendix: Yih’s solution 


Yih (1959) has considered the flow in a channel —1 < y < 1, x < 0, in which 
there is a sink at the origin and 


u—>cosinmy for «> 


so that = 
He gives the solution 


exp {(n?— })3 az} sin 


where A,(x,y) = it 
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We sum a series whose terms are asymptotically equal to A,. For n > «, 


mx 1 (Fes (-—1)” 


: 1 
A, (x,y) = sin nay + owns) 


n 8 
With = Jt, 
(+t) 
we have =—log(1Ft) <1,¢+ +1), 
n=1 


and related series can be ‘summed’ by integrating both sides from 0 to t. Hence 


\n 8 n(n—1) 8 128 4 }n(n—1)(n—2) 

= —log (1 —t) — log (1—t) +8} 
+ + log (1—t) — 
(1 +t) + 23] 

= (x,y), say. 

The function y—(4/7*)g has derivatives with respect to x and y up to the 
second order, for all x < 0 and all y, since the terms of its expansion are O(t"n-*) 
for n -> 00, with |t] < 1. Hence any singularities of / leading to infinite values of 
y or of its first and second partial derivatives must also be present in (4/77?) g. 
Near x = 0, y = —1, with 

zti=2*, 
the local singularity of y is therefore given by the (1+1t)*log(1+¢) term of g. 
Since z* = 0 is a stagnation point, we have 


= log 2* + O(|z*|?). 


This agrees with (11), for w = $7 on the streamline in question. 
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The reactant concentration spectrum in turbulent 
mixing with a first-order reaction 


By STANLEY CORRSIN 
Department of Mechanics, The Johns Hopkins University 


(Received 14 January 1961) 


The power spectrum of a passive scalar contaminant undergoing a first-order 
chemical reaction and isotropic turbulent mixing is deduced for three different 
spectral ranges: (i) the inertial-convective range; (ii) the viscous-convective 
and viscous-diffusive ranges for very large Schmidt number; (iii) the inertial- 
diffusive range for very small Schmidt number. The analysis is restricted to 
stationary, locally isotropic fields, and to systems so dilute that the heat of 
reaction has no effect on the reaction rate. 


Introduction 

A clarification and extension of the stationary, locally isotropic, turbulent 
mixing problem has recently been given by Batchelor (1959), and Batchelor, 
Howells & Townsend (1959). After preserving the Obukhov (1949)—Corrsin (1951) 
spectral result, G(k) ~ (1) 


with proportionality constant of order one, for the inertial-convective wave- 
number range, Batchelor focused attention on fluid+contaminant fields with 
very large ‘Schmidt number’ or ‘Prandtl number’, v/y7 > 1. vis the kinematic 
viscosity, 7 is the diffusivity, k is the wave-number, ¢ is the turbulent energy- 
dissipation rate, €, is the corresponding scalar contaminant field property, the 
rate of destruction of @ by diffusion. 

He pointed out that for v/y > 1 and large Reynolds number, there will be a 
viscous-convective wave-number range, 


in which the velocity spectrum is strongly affected by viscous forces, while the 
scalar field spectral transfer is predominantly convective. His analysis for 


k > (e/v3)t yielded yh] k\2 
G(k) = (3) 
B 


where k, = (e/v7?)t. This applies from the viscous-convective range (where the 
exponential is of order unity) all the way to k > oo. 

Then Batchelor, Howells & Townsend considered the opposite extreme of a 
fluid-contaminant combination for which v/7 < 1, in which there can exist an 


inertial-diffusive spectra range 
(i) 
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Here the turbulence spectrum is of Kolmogorov type 

F(k) ~ (5) 
with proportionality constant of order one, while the contaminant spectrum is 
strongly diffusive. For this range they deduced essentially that 


Lege! 


G(k) ke, (6) 

The present study is aimed at generalizing equations (1), (3) and (6) to include 
the effect of a first-order chemical reaction, assuming a system so dilute that 
associated temperature changes have negligible effects on the reaction rate and 
flow. Some preliminary results on decaying, reacting fields have been discussed 
(Corrsin 1958), and the formal consequences of the quasi-Gaussian hypothesis 
have been analysed (O’Brien 1960). 

The differential equation for concentration of a dilute contaminant under- 
going a first-order chemical reaction along with the convection and molecular 


diffusion is 
where C is a constant. With less dilute reactant, it will be necessary to include 
the heat of reaction and the temperature dependence of C. u is the turbulent 
velocity field, 7 is the (constant) diffusivity, ¢ is time, x is the co-ordinate vector. 
As pointed out earlier (Corrsin 1958; O’Brien 1960), this particular reaction-rate 
term is not spectrally selective, so in a freely decaying field it merely multiplies 
the mixing wave-number spectrum by an additional (exponential) decay in time. 
For stationary states, however, such simple arguments cannot be applied; 
we must pursue spectral elements in time as they migrate to ever larger wave 
numbers in the stationary spectrum. For the inertial range analysis of the 
turbulence itself this is essentially the approach used rather informally by 
Onsager (1949), arriving independently at Kolmogorov’s result F ~ k-3. For 
the viscous-convective spectral range in scalar mixing, this method of reasoning 
was introduced more rigorously by Batchelor (1959). Direct dimensional 
reasoning will not yield the explicit result of equation (23) because the reaction 
introduces an additional dimensional constant without increasing the number 
of dimensions. 
This paper is restricted to local isotropy for both velocity and reactant fields. 


Pure mixing in the inertial-convective range 

As a prelude to the following section, it is instructive to re-derive the k-beha- 
viour of (1) by a method like Onsager’s. We restrict our consideration to a k-range 
in which the transfer of @? spectral elements is purely convective. This means that 
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Employing the crude concept of a spectral ‘energy’ cascade in which the @? 
spectral elements ‘jump’ to successively larger k, we introduce the following 
assumptions. 

(a) The cascade wave-number sequence is a geometrical progression, so that 


Ak x k. (10) 


For turbulence dynamics this is suggested by the non-linearity of the Navier— 
Stokes equations; here we appeal to the quasi-non-linear traits of (7), due to the 
convective ‘forcing function’ coefficient u. 

(b) The part of @? transferred in a jump at wave-number k is G(k) Ak. We take 
it to be x G(k)k. 

(c) The characteristic time for each jump depends, of course, on the velocity 
spectrum F(k), and the simplest dimensional possibility is 


= k-3F-3, (11) 


The use of this form for the spectrally local convection time restricts us to the 
inertial part of the turbulence spectrum. 

In the inertial-convective range, the rate of flux of /? spectral content through 
any k must be independent of &, i.e. 


12 

dk\ 7) (12) 

Using (11), GF?k? = const. 13 
g ( (13) 


Using F ~ k-! as appropriate here, we find from (13) that 
G = Dk-, (14) 


with D a constant, consistent with (1). This approach does not immediately give 
the dimensional coefficients as in (1) and (5). They can be estimated by other 
means. For example, we may use the expression for diffusive destruction of 


spectral content, 7” 
= 2a dk. (15) 
0 


Assume that the order of the integral is given by the part between a wave- 
number ,k; characterizing the large structure and the wave-number (e/7%)t 
(for the cases v/y < 1). Then 


2yD| ks dk. (16) 

Integrating and neglecting oki. relative to (¢/y)}, we get 

D (17) 


consistent with (1). 


Reactive mixing in the inertial-convective range 
In the non-diffusive spectral range, (7) can be approximated by 


—CO,. (18) 
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6, represents any ‘narrow-band’ 0-field in the convective spectral range. Multi- 
plying by @,, and averaging, we lose the convective term by homogeneity and 
are left with the simple decay equation 


—2C06%.. (19) 
Integrating, we obtain the result 

On, = er. (20) 


The time rate of increase of wave-number for this reactively decaying spectral 
packet is still given by (10) and (11), since the first-order reaction has no effect 
on the geometry of the isoplethic surfaces, so 


(21) 


With (20), we see that (12) generalizes to 


d (Gk\ Gk dt d 


Using (21) for dk/dt, and (5) for F(x), this integrates to 


G(k) = Bi-Fexp | (23) 


which is shown qualitatively in figure 1. Evidently, for wave-numbers well 
above k, = Ce-}, the effect of reaction on spectral shape is negligible. &, locates 
the spectral region in which the rate of reactive loss is of the same order as the 
rate of convective spectral transfer. 

The dimensional coefficient B in this generalization of (1) or (14) cannot be 


estimated by a condition as simple as (15) because 20 | k?G'dk no longer repre- 
0 


sents the total rate of destruction of 62. Of far narrower applicability than (15) 
would be estimation of D in (14) or B in (23) by 


k* re) 
= | Gdk | (24) 
0 ok 


Here ,k, is the ‘energy-bearing’ wave-number; ,k* is the end of the inertial- 
convective range. The limitation is, of course, that the main contribution to this 
integral comes from the neighbourhood of ,/;, where G is unlikely to be well 
approximated by (14) or (23). The same limitation arises if we replace (15) by 
an expression for the total rate of destruction of 6?: 


= 29 | k2Gdk +20 | Gdk (25) 


ok* or 
| k2Qdk | Gdk. (26) 
oky oky 
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Having recorded these extreme doubts, but having no clearly better choice, 
we use (24) to estimate B. (26) is more tedious and not necessarily more reliable. 
Thus we write 
x B | k-bexp {3Cetk-B dk. (27) 
By integration, we get the final result 
— lexp {3Ce-5 — exp {3Ce-s 3}]. (28) 
Since ,k* > ,k,, we neglect the first exponential, getting the estimate 


~ exp { — 306-1 4kz3}. (29) 


log G 


FicureE 1. Qualitative sketch of reactant spectrum in inertial-convective range. 


To gain some appreciation of (28) or (29), consider its limiting form for C + 0. 
This is obtained most easily by putting (14) into (24), or by expanding (28) as 


a series, to give okt (20) 
Equating this estimate to the better one [equation (17)], we find that 


A similar approach to the inertial turbulent-energy spectrum gives the well- 
known estimate for dissipation rate in terms of large-scale parameters, 


ex (w)ik,, (32) 
omitting a factor of order unity. Putting (32) into (31), we get a form looking 
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If the large structures of velocity and contaminant fields are about equal, 
= O(u2)! ky. (34) 


In any case, we see that (29) represents the same kind of approximation as one 
used in obtaining (32). 


Reactive mixing in the viscous-convective range and beyond 


In the limiting case of mixing at very large Schmidt number, v/7, Batchelor 
(1959) pointed out the possibility of a locally isotropic viscous-convective range 
as outlined in the Introduction. His initial arguments are unaffected by the 
addition of our first-order, passive chemical reaction. 

This spectral range is characterized by the property that the ‘smallest eddies’ 
of the turbulent motion are much larger than the characteristic lengths of the 
sealar field, (¢/v3)! < k. Therefore, if we follow a small region of fluid, the turbulent 
straining action can be visualized as locally uniform for a limited time. This 
permits analysis of the history of a Fourier element of concentration field under 
the joint action of turbulent straining and molecular diffusion. Inference of 
general spectral behaviour from a single spectral point is permissible here 
(unlike analysis of the turbulence field) because the diffusion equation is linear. 

Batchelor chooses Cartesian axes which translate with a fluid element and 
are always alined with the local principal strain-rate axes. Calling these local 
principal strain-rates a, 2, y, meee (7) can be written as 

Presumably , / and y can be taken constant over distances comparable with 
the Kolmogorov microscale (v3/e)t, and over times comparable with the Kolmo- 
gorov time (v/¢)?. The subscript on 0, reminds us that we shall analyse the history 
of a single Fourier element of the full 0-field. 

Following Batchelor’s procedure, we get instead of his equation (4.6), the 
equation 


6,(x,t) = Ayexp sin (kz), (36) 
where ret, = (37) 


and A, is initial amplitude (see equation (4.4) in Batchelor’s paper). If we call 
0,,(x, t) the field with C = 0, and compare with Batchelor’s form for @,,, then 


0,({x,t) = 0,,(x,t)e-—. (38) 


To use such determinate results to infer the behaviour of a random 6-field 
randomly strained by isotropic turbulence, we argue that (36) applies to any 
Fourier element whose k > (e/v3)t . The magnitude of the maximum negative 
strain-rate is taken to be the root-mean-square value of strain-rate for the 
turbulence, i.e. we set vend (e/v). (39) 


Supposing that the 6-field is maintained statistically steady by continuous 
feeding into the large structure, we use (37) to estimate the kinematic history of 
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any narrow slab of the broadband 6-spectrum as it migrates to larger and 
larger k under the convective straining of the turbulence. Then (36) tells us how 
fast this slab loses @?-content during any time interval of spectral travel. 

The time required for a narrow slab to travel from k’ to k follows from (57). 
It is 


k 
given by Batchelor. The change in slab width as it goes from k’ to k is found to 
= Ok’ |k’. (41) 


During this time interval, the contribution to 6? carried by this slab decreases 
in the ratio given by the square of (36), with (40) for the travel time, or 


vA, = exp (42) 
To compare the actual spectral levels at k’ and k, we start with 
G(k) dk = ok’ (43) 
By use of (41) and (42), this leads to 


By inspection [or by using (44) we set up a difference equation, thence a differ- 
ential equation which can be solved], we arrive at the spectral form 


G(k) = exp (45) 


which reduces to Batchelor’s mixing form for C = 0. N = const. Introducing 
(39) for y, we derive the final result 


= Nk {— 2k (46) 


where k,, = (¢/vy?)!. The wave-number k, = C2e-} does not play a direct role in 
the shape of this part of the spectrum. It does, however, affect the entire local 
level through NV. 

In the viscous-convective range, k* < k < kp, (46) simplifies to a power law, 


G(k) x Nk-aracebe (47) 


We note parenthetically that (47) (and Batchelor’s limiting form for C = 0) 
can also be obtained by application of (22) with the spectral transfer time taken 
as T = const. = (v/e)}. 

Qualitatively, the combination of (23) and (45) gives a @-spectrum like 
figure 2. 

To evaluate V, we require that (46) [hence (47)| match the lower wave-number 
form, equation (23), at k = k* = (e/v?)t. This gives 


N = exp (48) 
or N = exp {30 (v/e)4}. (48a) 
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If we take (29) as a rough approximation for B, and neglect oki. relative to k*5 
in the exponent, (48) becomes . 
J 
It should be noted that the argument of the (constant) exponential can be 
written as {—3(k,/k,)3} [see remarks after equation (23)]. The ubiquity of the 
dimensionless reaction rate C(v/e)? is notable here, as in earlier results. 


C=0 ~k-3 
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FicureE 2. Qualitative sketch of reactant spectrum for v/7> 1. 


The relation (46) reduces immediately to Batchelor’s mixing case when we put 
C = 0. The corresponding simplification of the constant coefficient N follows 
by noting that (48) then reduces to 


N > Bk*-i, (50) 
while B>D (51) 
in accordance with (1). From (50) and (51), 

N +e,(v/e)3, (52) 


in essential agreement with the C = 0 theory. 


Reactive mixing in the inertial-diffusive range 

For v/n < 1, there exists an inertial-diffusive range in the 0-spectrum, its wave- 
number characterized by (4). For this case we can simply follow the procedure 
of Batchelor et al. (1959), with the addition of the reactive term as in equation (7). 
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The Fourier transform of (7) can be written as 


The simplifying theoretical arguments presented by Batchelor et al. (1959) 
can be summarized and adapted as follows. 

(a) Since the 0-spectrum doubtless falls off much more rapidly ie the 

u-spectrum in this range, the principal contribution to the integral must come 
at values of k’ of the order of and less than (¢/7%)!. It follows that over most of 
the spectral range of this analysis ki <k (54) 


in the important part of the integral. 

(b) The diffusion and reaction time for a Fourier component of @ around wave- 
number k is of order (yk? +C)-1, but the convective (quasi-source) effect has the 
inertial range characteristic time (ek?)-> [which is the same as equation (11)]. 
With k > (e/y)', the latter time is much the larger. Therefore the time-derivative 
term can be neglected in equation (53) when it is applied to the wave-number 
range of interest here. 

Neglect ¢P/ét in (53), then multiplying the remainder by its complex conjugate, 
and continuing to follow the Batchelor, Howells & Townsend arguments by 
using (54) to justify neglecting the statistical connexion between the (spectrally) 
widely separated factors in this integrand ({k—k’| > k’), we finally get the 
following extension of their equation (8): 


(nk? + C)PG(k) x F(k)/y, (55) 


where ¢# is the rate of destruction of @ by diffusion. Introducing (5) for F, 


we alrive at 
1 1 


| 56 
3 ki +07’ | (56) 


| 


which reduces properly for C = 0. With C + 0, it approximates the pure-mixing 
form in the wave-number range k > (C/7)!, which is a measure of the spectral 
location where reactive loss rate and diffusive loss rate are comparable. 

Since (56) has no arbitrary constant coefficient (except of order unity as a 
consequence of the various approximations en route), the matching of (56) with 
(23) at the wave-number ,k* = (e/y?)! just gives a relation among a number of 
characteristic constants in the problem. The matching leads directly to 


= EGE et 57 


With (29) as an estimate for B, this can be written as 


2 
1+0 exp{-se, | -1]}, (58) 
€ 


This approximation is not correct for C > 0. The non-reacting limit requires 
that the entire expression (28) be used for B. 
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The qualitative characteristics of reacting spectra for v/7 < 1 are sketched in 
figure 3. We notice in (58) that for this case the appropriate dimensionless reac- 
tion rate is C(7/¢)! instead of C(v/e)?. 
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FicurE 3. Qualitative sketch of reactant spectrum for v/7 <1. | 
Since (58) gives an estimate for the rate of diffusive destruction of #2, the total | 
rate at which (? is fed into the spectrum is just 
Gy = +2002. (59) 
This work was supported by Fluid Dynamics Branch, U.S. Office of Naval 
Research. 
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Theory of unsteady flow about thin cylinders in 
fluids of high electrical conductivity 


By LEON E. RING 
Cornell University, Ithaca, N.Y. 


(Received 10 June 1960 and in revised form 10 March 1961) 


A theory is developed for the incompressible flow of a fluid with high electrical 
conductivity about thin cylinders (airfoils) in non-uniform motion. A uniform 
magnetic field is applied parallel to the free stream and solutions are obtained 
subject to the restriction of small perturbations. The effects of viscosity are 
included, for the most part, only through the use of the Kutta condition, where 
applicable, for lifting airfoils. The validity and range of applicability of the 
infinite-conductivity assumption are determined on the basis of an order-of- 
magnitude analysis; the general character of the flow is discussed at length. 

The flow-field for infinite conductivity is changed from the non-magnetic 
case only through the new transport speed of vorticity; the forces on the airfoil 
are changed due to surface currents. For the case of the Alfvén speed less than the 
free-stream speed, the airfoil lift and pitching moment are given in integral 
form for general unsteady-airfoil motion and are given in closed form for har- 
monic ocsillations. The forces at moderate frequencies may be larger than in the 
corresponding non-magnetic case. The response to a unit-step change in the 
downwash is studied and the asymptotic form of the lift is obtained for small 
and large time. 

For the case of the Alfvén speed greater than the free-stream speed, vorticity 
and current are shed from both the leading and trailing edges. Therefore the 
extension of the usual Kutta condition is not obvious. It is shown that if finite 
viscosity and/or conductivity tend to remove the trailing-edge singularity, the 
flow is unstable and no steady flow can be obtained. 


1. Introduction 

The steady magnetohydrodynamic flow of an incompressible fluid of large 
electrical conductivity about thin cylinders (airfoils) has been studied by Sears 
& Resler (1959). In the limit of infinite conductivity with a uniform magnetic 
field applied parallel to the free stream, the flow field is found to be irrotational 
and current-free and is in fact identical with the corresponding non-magnetic 
flow. The magnetic field inside the body is identically zero, so that the pressure 
on the body surface is altered due to the implied current sheet on the body surface. 
For large but finite conductivity, this current sheet diffuses into a thin layer. This 
inviscid ‘magnetic-boundary-layer’ has been studied by Lewellen (1959) and 
is found to have a thickness of the order R;,? , where R,,, is the magnetic Reynolds 


number. 
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The results of these investigations are of sufficient general interest to suggest 
an analogous study of unsteady flow. In particular, since thin cylinders in 
asymmetric flow generally produce circulation and lift, one is led to consider the 
effects of unsteady circulation and vortex shedding. Thus, although ostensibly 
concerned with ‘airfoils’, the present paper is intended to cast light on phenomena 
that can be expected to appear in a rather general class of unsteady magneto- 
hydrodynamic flows around bodies. The term ‘airfoil’ is used since it suggests 
the thin cylinders being considered and since the methods of classical thin-airfoil 
theory are applicable to the present problem. 

A free-vortex element in an incompressible and infinitely conducting fluid will 
travel at the Alfven speed relative to the flow. For the aligned-fields case, this 
says that shed vorticity will travel parallel to the free stream, but will have a new 
transport speed. This suggests the use of the techniques of classical unsteady- 
airfoil theory, as given by Karman & Sears (1938), with modifications for the new 
vortex-transport speeds and the airfoil forces due to surface currents. This is 
done in the present paper, so that the Sears & Resler (1959) steady solution with 
aligned fields is extended to include unsteady airfoil motion. Of course, this is 
based on the assumption, made by Sears & Resler and others, that the correct 
steady-flow pattern has undisturbed conditions far from the body. Other 
hypotheses, leading to grossly different flow patterns, have been suggested by 
Stewartson (1960). 


2. Basic equations 

The equations which govern the motion of an incompressible fluid with finite 
viscosity and scalar conductivity have been given by a number of authors 
(e.g. Resler & Sears 1958) and can be reduced to the following dimensionless form: 


Dq 1 22) — +2 7 72, 
V.q=0, (3) 
V.B=0. (4) 


The equations have been non-dimensionalized using the free-stream speed U*, 
the applied magnetic field of induction Bj, and the body length Lx. M.x.s 
rationalized units are used, with an asterisk denoting dimensional quantities. 
The three parameters are defined by 

U*L* 

R,, = a= By 


R,, is the magnetic Reynolds number and R, the viscous Reynolds number. 
The parameter «? is the ratio of the magnetic pressure to the dynamic pressure 
in the free stream; it is also the square of the ratio of the speed of an Alfvén-wave 
in an infinitely conducting medium to the free-stream speed. The fluid is assumed 
to have uniform permeability ~* and scalar conductivity o*. We will also use the 
vorticity 2 = curl q and the current density j = curl B. 
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Now consider the two-dimensional flow where small perturbations are super- 
imposed on a uniform free-stream and a parallel magnetic field. This is repre- 
sented by q = (1+u,,u,,9), u,,u, <1; 
B = (1+6,,b,,0), 6,,b,<1. 


For the most part, it is assumed that the flow field is described by taking zero 
viscosity and infinite conductivity in the linearized form of equations (1) and 


(2), which become 
dp +a% )= (5) 
da’ 


cb cb ou 
ct 


(6) 


3. General character of the flow 
Steady flow 


Before seeking a solution, it is well to establish the general character of the flow 
on the basis of order-of-magnitude arguments. In particular, the validity and 
range of application of the infinite-conductivity assumption (to be used later) 
will be established. First, let us consider flow with R,, < R,. 

By comparing the first and last terms of equation (2), one may write a charac- 
teristic magnetic diffusion time as 7,, = R,,65,, where 6,, is a characteristic dif- 
fusion distance or, in the present case, the current-penetration depth normal to 
the body surface. In steady flow, the time available for diffusion is the time 
required for current to travel the length of the body. Since the current travelling 
with the upstream Alfvén mechanism remains in the neighbourhood of the body 
longer and thus allows greater time for diffusion, the characteristic convection 
time is taken as 7* = L¢/|U%—Bgu*-p*-3|, which in dimensionless form 
becomes 7, = 1/|1—a|. Equating these two times, the penetration depth for a 
convection-diffusion balance is given by 6,,, = (R,,|1—a|)~*. This agrees with 
the result of Sears & Resler (1959), namely that for infinite conductivity the 
flow field is irrotational and current-free. It also agrees with Lary (1960), who 
found the penetration depth to be |k|-} = {4R,,|1—«2|}-4, since }(1+«) is of 
order unity. 

The strength of the magnetic field inside the body can be estimated by noting 
that the magnetic flux inside the body must be equal to that which has diffused 
from the fluid. The flux which has diffused from the fluid can be approximated by 
6,,\B, where AB is the change in the field strength across the diffusion region. 
The flux inside the body can be approximated by eB,, where ¢ is the body thick- 
ness ratio and B, is a measure of the field strength inside the body. Equating 
these, B, is given by B,= 6,,AB/e. If 6, < ¢, i.e. if the diffusion depth is much 
less than the body thickness, then B, must be nearly zero since AB can be no 
larger than order unity. This gives an effective ‘no-slip’ condition on B and 
results in a ‘magnetic-boundary-layer’, analogous to the usual viscous boundary 
layer. If 6, > €, By can be expected to be of the order of the free-stream value 
and AB a small quantity. This agrees with the Sears & Resler result that 
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the magnetic field is zero inside the body for infinite conductivity, and with 
Lary’s result that the magnetic field is only slightly perturbed inside the 
body when ¢ < {}R,,|1—22|}-#. 


Unsteady flow 


Now let us extend the above results to unsteady flow. Consider the pitching 
motion of athin cylinder (insulator) which is described by the circular frequency w. 
The unsteady characteristic time is taken as the ‘period’ 7,, = w—1. The pene- 
tration depth for a balance between diffusion and an unsteady input is found 
again by equating times to give 4,,,, = w}R,,?. 

When 6,,,, > 4,,,, the unsteady current component is carried away by the 
Alfvén mechanism before it can diffuse to the depth 4,,,,,, so the essential balance 
remains between convection and diffusion. With this balance, the general flow 
character (except for the wake) can be expected to be the same as for steady 
flow. 

For 6, < On. the unsteady current component is restricted to a smaller region 
than the steady component. Although there is a small-scale balance between 
diffusion and unsteady input, the large-scale balance remains between convec- 
tion and diffusion. There are three cases depending upon the body thickness- 
ratio ¢. (a) When 0,,,, < 6,,,<eé, both mechanisms tend toward a magnetic 
boundary layer. (b) When 4,,,, < € < 4,,,, the unsteady current gives a magnetic 
boundary layer,while the steady solution gives a high degree of diffusion. This 
could be achieved for moderate values of w with a high R,, and « sufficiently 
close to unity. This flow can be described by saying that the magnetic field is 
effectively ‘frozen in the flow’ during an oscillation period but can diffuse a 
large distance during the time that the current stays in the neighbourhood of 
the body. (It is also ‘frozen out of the flow’ or ‘in the body’ during an oscilla- 
tion period.) (c) When e< 4,,,,, < 6,,,, both mechanisms allow a high degree of 
diffusion. 

For the slender bodies being considered, the boundary conditions at the body- 
fluid interface can be satisfied by discontinuities across a singularity sheet, 
so the solutions are linear and superposition of solutions is possible (for thickness, 
camber, unsteadiness, etc.). That this is true for the limiting cases discussed 
above follows directly from the work of Sears & Resler (1959) and Lary (1960) 
except for the special case 6,,,,, < € < 6,,,. That it is also true in this case can be 
seen by noting that the steady solution has only small perturbations everywhere. 
Therefore the unsteady solution, within the linear theory, sees only an undis- 
turbed flow. The boundary conditions for the solution of the unsteady flow field 
can thus be satisfied by a singularity sheet. The required unsteady current-sheet 
strength corresponds to the jump from the unsteady flow to the steady magnetic 
field inside the body, since the field strength inside the body cannot change 
during an oscillation period. 

From the above, it is now possible to answer the question of the applicability 
of the infinite-conductivity theory. For the calculation of thickness effects, in 
general it is required that the magnetic field be approximately zero inside the 
body; this requires that ¢ > (R,,|1—a|)-4. However, the unsteady infinite- 
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conductivity lifting theory only requires that the unsteady current be confined 
to a thin layer; thus infinite-conductivity theory is valid when either 


1<(R,,|1—a|)* or 1<(oR,,)}. 


A case of particular interest occurs when 1 < (wR,,)! and 1 > (R,,|1—a|)}. 
The solution for this flow is given by the supposition of a steady contribution, 
calculated on the basis of Lary’s (1960) theory, and an unsteady contribution 
calculated on the basis of the infinite-conductivity theory of the present paper. 


4. Conservation laws and Kutta condition 

Conservation of vorticity 
Consider a body at rest which starts to move at time ¢ = 0. At any finite time f, 
all current and vorticity must be in a finite region of the plane which can be 
enclosed by a curve C, where C is convected with the local flow. Since Q and 
j are assumed to be zero along C, Kelvin’s theorem can be applied to the curve C 
and one can immediately write that 


= (7) 


where the integral is taken over all the vorticity, that is, the net vorticity in the 
system remains zero. 


Conservation of current 


A similar constraint will now be obtained for the current. Consider the contour 
C to be a circle of radius r about the body, and integrate eyuation (2) around C. 
This gives P 
B.al =| (8) 
dt Cc Cc Rn 
Since there is no current on C’, V?B = 0 on C. The net vorticity was shown to be 
zero, So q = O(r-?) on C. Therefore the right-hand side of (8) is O(r-!), which can 
be neglected for a sufficiently large contour. Then, applying Stokes’s theorem 
and integrating over time, we have that 


| jdady = const. = 0, (9) 


where the integral is taken over all current, and the integration constant is 
evaluated at a time before the start of the body motion. This says that the net 
current in the system remains zero. 


+ One must be a bit careful about what is meant by a two-dimensional problem, since 
all currents must close somehow. If there is a net current flowing, there is some question 
if the closing of the current (in a large wire loop, for instance) can be neglected. In the 
present work, the net current is zero. Therefore the current elements can close through 
conducting side walls, which would give an effect analogous to tip effects of classical-airfoil 
theory. 
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Kutta condition 

In classical two-dimensional thin-airfoil theory, the undetermined cyclic constant 
(circulation) is specified, to make the solution unique, by using the Kutta con- 
dition. The criterion for small, but finite, viscosity has been discussed by a 
number of authors, including a recent comprehensive discussion by Sears 
(1956). It is because of the fact that vorticity convects with the flow (and diffuses 
relative thereto) that the boundary-layer vorticity leaves the airfoil at the trailing 
edge; therefore the Kutta condition must be applied at this edge. 

For a magnetic boundary layer with «? < 1, all vorticity and current must 
leave the airfoil from the trailing edge. Thus the Kutta condition must still 
be applied at the trailing edge. For « > 1, the extension is not clear. It can be 
shown that if the angle of attack is suddenly changed, there is an immediate 
flux of vorticity and current from both the leading and trailing edges. Therefore 
it seems that, at least in the unsteady case, the criterion should include the 
leading as well as the trailing edge. Since we lack such a criterion, which must 
come from more detailed boundary-layer considerations, this case will be treated 
without applying a specific condition. 


5. General unsteady airfoil motion 
Fluid model 


Consider the unsteady two-dimensional flow about a thin cylinder (airfoil) 
for arbitrary values of the parameter a?. The airfoil is taken to have zero thick- 
ness and to lie near the x-axis in the range [0, 1]. If the fluid has sufficiently high 
electrical conductivity and low viscosity, all current and vorticity are confined 
to thin wakes and the wakes can be represented as singularity sheets. 

The fluid model used for the solution of the flow field is shown in figure 1, 
and consists of a current and vortex distribution in a downstream wake, over the 
airfoil, and when applicable, in an upstream wake. The current and vorticity in 
the wakes can be related by the solution of equation (5) and (6) with the con- 
dition that the wakes are force-free. This gives the well known Alfvén-wave 
relations, eM(a,t) = ae? (a,t) = 

t) = = fyft—a/(1 
where € is wake vorticity, e® is wake current, and the subscripts 1 and 2 refer 
to upstream and downstream waves respectively. These have velocities and direc- 


tion cf travel as shown in figure 1. For any finite time, the wakes will be of finite 
length and all current and vorticity will lie in the interval —R < # < R. 


(10) 


Potential relations 


The flow outside the wakes can be given by a perturbation magnetic potential 
¢” and a perturbation velocity potential 6%. A relation between the potentials 
is obtained by substitution into equation (6). Integration of the «-component 
with respect to a yields the result 
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where the integration constant is found to be zero, by evaluation at x = x. é 
is now given explicitly in terms of 6 by integration of the above result, following 
a fluid particle. This yields the integral 


t AAD (> — 
t) = | dg (ita) 


=] 
a. 
° 


Velocity (1 +2)U, 
Uo Airfoil 
£ € ® 


EDEDEDE * 


Vortex _/ Velocity (1-a)U, 


Current 
(a) 
U, 
_* End of wake Airfoil End of wake 
-R ve, t) R 
Velocity Velocity +1)U, 
Vortex 
Current 


(b) 


FicureE 1. Fluid model for a lifting airfoil which has been in motion a finite time, with 
(a) a? < 1, (b) a > 1. Wake elements have been displaced from x-axis for clarity. 


The potentials are taken to be identically zero at tf = — oo. The magnetic and 
velocity potentials must satisfy equation (11) at every point in the flow field. 
Therefore the discontinuities must also satisfy this relation, and the current 
distribution y® can be expressed in terms of the vortex distribution y® in the 
same form as equation (11) (with no y dependence). The vorticity and current 
in the wake are expressed as 


4 (1), (2) x,t => (2) = (2) + 6M), (2) x,t 12 
Y 1 2 


where the Alfvén-wave relations are given in (10). 
The circulation and total surface current are related to the distributions by 


1 
TO. = | ys (ar, t) da, (13) 
0 


where the circulation is positive in the counter-clockwise sense. The total surface 
current can be expressed in terms of the vortex distribution by using the current 
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equivalent of equation (11) in equation (13) and the functional form of the wake 
from equation (10). Thus we derive the result 


T(t) yng, ef) (0-, 9) dy. (14) 


The first term on the right-hand side is the vortex distribution integrated in a 
time-lagging manner; the second is proportional to the vorticity shed from the 
leading edge in the previous unit of time. These indicate the general time-lagging 
character of the magnetic field to the velocity field. 


The analogous Wagner integral equations 


In classical unsteady-airfoil theory, the airfoil motion is related to the wake 
vorticity through the so-called Wagner integral equation. If we follow the tech- 
niques of the classical theory (e.g. von Karman & Sears 1938), we obtain the 
generalization: 


=f a. (15) 


wake 


Here the quasi-steady circulation is given by 


= 2 (4) dé, (16) 


and the airfoil surface is given by Y (x,t +a" is the linearized convective 


derivative. ['{) represents an arbitrary circulation; if the Kutta condition is 
applied at the trailing edge, then ['y” is identically zero. The integration over the 
wake includes the forward as well as the rearward wake. It should be noted that 
the conservation of vorticity, equation (7), is used in the derivation of (15). 

The total surface current is expressed in terms of the vortex distribution by 
(14). Using this, an analysis similar to that leading to equation (15) yields the 


result 
1/ é dx 

l-a 


nan, 07) 
t—-1 


where the quasi-steady surface current is given by 


7 x Dt 1—£} 
and the surface current corresponding to the arbitrary circulation is 
dx 
PO | (19) 


Equation (17) is the analogue of the Wagner integral equation for current. 
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The problem is now reduced to the solution of the integral equations (15) and 
(17), with the quasi-steady terms defined by (16) and (18), with the functional 
form of the wake given in (10) and given some law to specify the circulation I'\)) 
(e.g. the Kutta condition). 


The lift 


Consider the lift due to the vortex and current distribution over the airfoil and 
in the wake. Since the wake is force-free, it can introduce no net force. Divide 
the lift into that due to the vorticity and that due to the Lorentz force. The lift 
due to the Lorentz force is 


AL* = B* | | jtdx* dy*, (20) 


where the integral includes all current. This integration jusi gives the total cur- 
rent in the system, and, by equation (9), this is identically zero. Therefore the 
lift on the airfoil is equal to the force on the total vorticity in the system. The 
force on the vorticity can be obtained by calculating the impulse of the vortex 
system, as given by Karman & Sears (1938); this gives 


f* d 


— —~1)y@ | 2 
CO) = = (21) 


-R 
The lift is treated following the method of Karman & Sears with modification 
for the change in vortex transport speeds, used by Ring (1960). The result, 
including arbitrary circulation and a possible forward wake, is 


(Ie 4) —2— 
€ (§, ¢) 1 


C,(t) = A, OW) | 


wake 


wake 1 


1 
Az = 24 (23) 
dt 0 
Here y{) is the quasi-steady vortex distribution over the airfoil, 
2 t){ & \* dé 


From equation (22), the lift can be compared to that of classical unsteady- 
airfoil theory. The apparent mass and the quasi-steady lift (the second term) 
are identical with classical theory. The third term represents the effect of an 
arbitrary circulation and is identically zero for a trailing-edge Kutta condition. 
The effect of vorticity in the wake (the fourth term) is formally the same, but is 
changed due to the new transport speed of vortices. The last term is new and 
might be called the magnetic term. The kernel of the magnetic term is seen to 
approach unity far from the airfoil, unlike the kernel of the wake-vorticity term 
which approaches £—!. This means that the magnetic term takes into account all 


where A, is the lift associated with apparent mass and is given by ar 
= 
_| 
8) 
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current inthe wake and gives a contribution even after the flow has become steady 
near the airfoil. In fact, it can easily be shown that the magnetic term includes 
the effects of surface currents in the steady solution and gives rise to the (1 — «?) 
factor of Sears & Resler (1959). 


The pitching moment 


As in the lift calculation, we divide the pitching moment into that due to the 
vorticity and that due to the currents. The pitching moment coefficient about the 
airfoil midchord (x = }), measured positive in the clockwise sense, due to the 
Lorentz force can be written as 


AM* “a 2 
The current can be expressed in terms of the vorticity in a manner similar to the 
potential relation (11a), namely 


R t ay 
AC,, 2a (x— 4) ~— (€+a-t, £)dé. (26) 

We differentiate this with respect to time and then integrate by parts with respect 
to x; then, by use of (11a), this becomes 


(27) 


The first two terms on the right-hand side are zero, since they represent vorticity 
and current beyond the ends of the wake; the last term is the net vorticity minus 
the net current in the system and therefore is zero. By integration over time, the 
incremental pitching moment is found to be given by 


AC, 


‘nn = constant = 0, (28) 


where the constant is evaluated at a time before the start of the airfoil motion. 
Therefore the surprising result is obtained that the pitching moment due to the 
total current distribution is zero. The lift due to the current is zero essentially 
as a consequence of the conservation law for current; however, it is not obvious 
why the pitching moment should have been zero. 

The pitching moment due to the vorticity can be found by using impulse 
techniques. Again following the method of Karman & Sears as modified by Ring, 
the pitching moment can be expressed as 
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where 4, is the pitching moment related to apparent-mass effect and is defined 

by ap 

By = al t) dx. (30) 
0 


The first three terms are recognized as those from classical unsteady-airfoil 
theory. The wake-vorticity term is formally the same but is changed because of 
the new transport speeds of vorticity. The last term is new and might be called 
the magnetic term. 

As the duration of the airfoil motion becomes large, and therefore the current 
is far from the airfoil, it appears that the magnetic term diverges, since the 
integrand approaches fe, t). However, that this is not actually the case can be 
seen by subtracting the magnetic moments due to surface current on the airfoil 
(expressed in terms of the wake from equations (25) and (28)). In fact, the 
magnetic term again gives rise to the (1 — «?) factor of Sears & Resler (1959). 


6. Harmonic oscillations, «? < | 
Calculation of wake distributions 
Consider an oscillating airfoil with the upwash on the airfoil given by 


DY (z,t) 
Dt 


= g(x)e™, (31) 


where g(x) is an arbitrary function. The frequency w is taken such that Ziw > 0; 
this corresponds to an airfoil oscillating with increasing amplitude. Clearly 
#iw can be put equal to zero in the final results. 

Following §4, the Kutta condition is applied at the airfoil trailing edge, since 
a2 < 1. This means that [' must be identically zero. 

The quasi-steady circulation and surface current can be expressed as 


= GU Tr?) G2 (32) 


where G® and G® are complex in general and are given by using equation (31) 
in equations (16) and (18). Thus we obtain the relations 


1 
GM = 2 g(é) (45) dé, (33) 


From the functional form of the wake given in (10), the wake distributions can 
be written as 
t) = t) = gy exp {iw[t — (w—})/(1—a)]}, 


(1) (2) (35) 
t) = — (a, t) = geexp {iolt— (e—3)/(1 + 


where g, and g, are complex constants to be determined. Using equations (32) and 
(35), equation (15) can be reduced to 


— 26 = KP + KP] KP + KP], (36) 
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where the integrals have been identified as Bessel functions; K,,(z) is the nth- I 
order modified Bessel function of the second kind. The superscripts are a short- 
hand notation to indicate the arguments, as shown below: 


1) w whi 
(37) | 
The factor (1—.«)/|1—a|, which is unity in the present case, is included to make 
the definitions useful for «2 > 1. Similarly, equation (17) can be written as The 
ane 
— 24° = g, — a), (38) mor 
where D is defined by T 
D(w, «) = = exp | ae 
(1—«) Cas 
(39) the 
The function D has been evaluated in Ring (1960) and is late 
D(w, = KP +a KY (4C) 
where J,,(z) is the nth-order Bessel function of the first kind. The implied argu- a 
ments of Jjand J, are $v. Thesimultaneous solution of equation (36) and (38) gives i 
as 
g, = Do, KP] (41) 
Dio, (KP + KP] — a) 
where G = G2/G%, (42) | 
The constant g, is given by taking « -> —« in the expression for g,. This essen- C 
tially represents the solution of the flow, since the lift, pressure distribution, etc., ? 
can easily be expressed in terms of G®-® and g, 9. 
The 
Force calculations al 
The lift has been given in equation (22); with the use of (35) and (41), this can be T 
Cr(t) = 2), (43) | 
where wit] 
(hic _ (l-a) KP{Do, — a) + — (1+ a) KP{D(o, «) — G[Ky + KP} the: 
Div, — 2) {KY + KP] Dio, a) [KP + KP] 
(44) tenc 
T (Siw, x) is the so-called Theodorsen function of classical unsteady-airfoil theory, 
modified to include the magnetic effects. (It should be noted that w is the Tak 
reduced frequency based on the full chord, as opposed to the usual practice of 
basing it on the semi-chord.) The most important change in the Theodorsen 
function is the appearance of G. Thus, while the classical Theodorsen function 
depends only on frequency, the modified Theodorsen function depends upon the 
and 
mode of oscillation (and on «) as well. 
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In a manner similar to the life calculation, the pitching-moment is found to be 


1 
C,,(t) = Ay +{ — 1) (a, t) dx + AT Q(t) (45) 
0 
where 
(1—a) — (4a/iw) KY?) [D(o, — -—G{KP + KP}] 


D(o, — a) [KY + KY] «) [KP + K®) 


The function S(37w, «) includes the wake and all magnetic effects. The appear- 
ance of G again shows that the function is not universal but depends upon the 
mode of oscillation considered. 

Three G’s are of particular interest, namely, 


Gy = lo, = + 24/0. (47) 


Case I corresponds to vertical airfoil oscillations, case II to a downwash linear 
about the mid-chord, case III to a downwash distribution linear about the quar- 
ter-chord. The first two are of interest because of their mathematical simplicity, 
the last because of a special high-frequency property which will become apparent 
later. 
Low-frequency limit 

The low-frequency expansions of 7’ and S are obtained by expanding @ in a 
power series and by using the small-argument Bessel-function expansions. 
After straightforward but lengthy algebraic manipulation, the force coefficients 
are found to be 


C,(t) = (1 +a?) Ay — (t) (1 — 2?) + 2”) ion fo + +2?) 
x {y+ hin — (1—2*)} + O( In? Jw)}, (48) 


1 
C,,(t) = Ay + (1—a2) | (20 — 1) (x,t) de — +2) ion Jo + 
0 
(49) 


The apparent-mass term is included in the pitching-moment although it is of 
order w and some terms of order w have been neglected. 

The first thing to note from these expressions is that as w — 0 the force co- 
efficients become (1 —«?) times the values for non-magnetic flow in agreement 
with the Sears—Resler result. The coefficient of the (wInw) term is increased by 
the factor (1 + «?); the coefficient of w is also increased because of magnetic effects. 
Thus, although the steady forces are reduced by a factor (1 — «?), there is a general 
tendency for the unsteady contribution to be increased by the magnetic effects. 


The limit a > 0 
Taking the limit « > 0 (i.e. the magnetic field approaching zero), it is seen that 
T (iw, a) = ), (50) 
and S(iw, a) = = + O(a?) (51) 
Ky (iw) + Ky (tw) 
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These limiting forms are identical with classical theory. The coefficient of the 
a* terms can be obtained; however, it turns out that the frequency dependence 
is of such a complicated nature as to yield little information. 


High-frequency limit 
Asw—> %, it is convenient to represent G by an asymptotic expansion of the form 


Using this expansion of G, the asymptotic expansion of 7’ is 
T (diw,a) ~ 2G) e-iv + {1 + 20? — 
+ ia®[1 — 5G — 46-]e- 4 — 2G) (53) 


This expansion is valid for —7 < argw < 0; the top sign is to be used for 
argiw > 0 and the bottom sign for arg iw < 0. The expansion of S is similar. For 
z= 0, T is asymptotic to the constant 3, as is known from classical theory. 
However, the general case represents a decidedly different behaviour. If viewed 
in the complex plane, as w becomes large (and real), 7’ describes a circle of radius 
1x2(1—2G) about the point 4. It approaches this circle in an oscillating manner, 
as seen from the second term in the expansion. Thus, while the steady-flow 
theories of Lary (1960) and of Sears & Resler (1959) predict small forces in the 
neighbourhood of « = 1, the present theory predicts large forces due to unsteady 
effects. 

Case IIT (equation (47)) has the special property that (1 — 2G) = 0; in fact, it is 
the only linear downwash distribution to have this property. The reason why this 
is true will be pointed out in the discussion of the Wagner problem (§7). 


The limit a—>1 


There is some question as to the validity of the solution for « = 1 (e.g., the Kutta 
condition). However, the limiting solution is of interest since it will shed some 
light on the solution for « near but not equal to unity. For « > 1-, we obtain 


J,-G 


= K) 
T (}tw,1-) = KY +i, (54) 
4 
S(dtw, 1-) = T(diw, 1- (55) 


The behaviour of these functions for large frequencies is similar to the general 
case, except that the radius of the asymptotic circle has been maximized. It is 
not valid to take the limit w — 0 in these functions, since for « near unity we 
required that 1 < (wR, ($3). However, for large R,,,, @ can be quite small. Thus 
we can consider small frequencies, keeping in mind the physical requirement that 
w must be positive. Using the low-frequency expansions of equations (54) and 
(55), the force coefficients are 


C,(t) = 24,+O0(v?), Cm/(t) = O(0). (56) 


Thus for low frequencies and « near unity, the forces on the airfoil are small. 
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The circulation for 1~ and small wis 
1 
Pat) = tio | (22 —1) (a, (57) 
0 


so that the circulation is small. For zero circulation, the flow pattern must 
correspond to classical hydrodynamic flow with a relaxed Kutta condition. There- 
fore the flow pattern approached by taking w > 0 with « = 1~ has zero lift with 
zero circulation, as opposed to the limit in the Sears—Resler theory (and in the 
present theory), « > 1, with w = 0, which gives zero lift with 1 = 0. Following 
§3, both limits may be valid for &,, = 0%, withneither valid for R,, < «2. However 


the former limit can be approached to any desired degree of accuracy by taking 
R,, arbitrarily large but finite. 


7. The Wagner problem, «? < 1 
Formulation of problem 


The classical problem of the response to an instantaneous change in the down- 
wash distribution on a thin airfoil is referred to as the ‘Wagner problem’. The 
downwash for this case may be represented by 


= g(x)H(t), (58) 


where H(t) is the unit-step function defined by 


0, t< 0, 


HW) = t>0 


The downwash of equation (58) can be obtained by Fourier superposition from 
equation (31). Similarly, all quantities from the harmonic-oscillation solution 
may be superimposed to give the corresponding quantities caused by the in- 
stantaneous change of downwash. 

The quasi-steady circulation is obtained from equation (32) as 


T(t) = GMH(t) (59) 


The lift, obtained from equation (43) is 


1 
O,(t) = d(t) | (22 —1) 0+) dx — 27 (t) W(t. x), (60) 
0 
est 
where W(t.a) = mil... T(28, a) —ds, (61) 


and 6(t) is the Dirac delta function. The integral has been expressed as a Laplace- 
transform by taking s = iw and 2, to lie to the right of any singularities. The 
apparent mass has given an impulse lift at ¢ = 0; the lift at later times is all due 
to the modified Wagner function W(t,«). Let us confine our attention to this 
function. 
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W(t, a) for small times 


The standard method of obtaining the behaviour of W for small times is to obtain 
the asymptotic expansion of 7 for s > 00, and then relate this to the asymptotic 
expansion of W for ¢t + 0. The expansion of 7'(}s, «) for s > o is given in equa- 
tion (53) and contains exponential terms. It has been shown by Ring (1960) 
that the e~* terms give a contribution like exp (— 1/t) for small times. Therefore, 
using known results (e.g. Doetsch 1950), the expansion of W for small times is 
given, from equations (61) and (53), as 

W(t,x) (62) 
The first term in W remains equal to 3, unaffected by the magnetic field. The first 
magnetic effects are proportional to time and, since G* has appeared, depend 
upon the downwash distribution considered. 

The result given in equation (62) will be valid even for small £,,, provided the 
time is taken so small that the current is still restricted to a thin layer. Following 
$3, this will be true for t < 10-?R,,, where the limit on the current-penetration 
depth has arbitrarily been taken as 6,, = 107}. 


The logarithmic singularity 
From the asymptotic expansion (53), it is seen that convergence of the integral 
in equation (61) is obtained for all time except at ¢ = 1 where e* balances the 
factor e~’ and gives rise to a logarithmic singularity. By extracting the coefficient 
of the singularity, W is found to be given by 


W(t,a) = —2G)In 


+0(1). (63) 


This represents a marked change from the classical Wagner function, which has 
no singularities. The singularity is seen to arise from the high-frequency oscil- 
lating behaviour of 7’; having connected these phenomena, it is now possible to 
examine their origin. 

The logarithmic singularity in W appears at time t = 1. At ¢ = 1, the fluid at 
the trailing edge is that which was at the leading edge when the downwash was 
changed. Thus we are led to the examination of the flow near the leading edge at 
high frequencies. A calculation of the vortex distribution yields 


t) = — (t) {(1 —x)/x}4 {(1 — 2G) + O(1/w) + O(x)}. (64) 


Thus the coefficient of the high-frequency oscillations of T (tiv, «), the coefficient 
of the logarithmic singularity of W(t, a), and the coefficient of the leading-edge 
vortex singularity at high frequencies all contain the factor (1—2G@). Thus it 
appears that the leading-edge vortex singularity interacts with the magnetic 
field to produce the observed behaviour of 7’ and W. 


W(t, «) for large times 
The asymptotic expression of 7’ valid for |arg s| < 7 is 


T(4s,a) ~ (65) 
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By setting the denominator equal to zero, 7'(}s,«) is found to have infinitely 
many simple poles to the left of, and bounded away from, the imaginary s— axis. 
It is shown by Ring (1960) that the integration in equation (61) along C; and 
C; (figure 2) is o0(1) as n > oo fort > 1. Therefore (61) can be written as 
1 


st st 
Wit,a)= x) ds Res _ (66) 


where the summation is to be made over all poles to the right of C,. The summation 
can be shown to converge for ¢t > 1, by using the asymptotic formulation at the 
poles. 


\ A A 
| 
x | 
| (2n+1l)a+in C; Complex 
A s-plane 
C; 
\ | 
\ 
/ | 
4 | 
4 
/ | 


FicureE 2. Integration contours in the evaluation of W(t, a). 


It is now assumed that 7'(3s,«) has no poles for #s > 0, i.e. that there are no 
unstable or neutrally stable solutions with ['f? = 0. By using the expansion of 
T for small s, the integral in equation (66) for large times can be given directly 
(e.g. Doetsch 1950). Hence, 

W(t,a) = (67) 
The summation term in equation (66) does not contribute under the asumption 
that there are no poles for Zs > 0 and because the effect of any poles with Zs < 0 
(which do in fact exist) is exponentially small. In this connexion, it should be 
recalled that these poles are bounded away from the imaginary axis. 

Equation (67) is seen to reduce to the classical case for a? = 0. As time goes 
to infinity, the solution approaches the Sears—Resler solution; the time to reach 
the steady solution is increased by the factor (1+). This effect might be ex- 
pected in view of the lagging of the magnetic field. 


2 
8. Flow with a* > 1 General relations 


The theory for general unsteady-airfoil motion has been formulated in §5 in 

terms of an arbitrary circulation I’), where ') is zero for a trailing-edge Kutta 

condition. For a? > 1, the Kutta condition is not directly applicable; also, 
28 Fluid Mech. 11 
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the proper extension is not obvious, as noted in §4. Without a law tospecify Ty), 
the flow cannot be determined uniquely. However, some general information 
regarding the stability of the flow can be obtained. For this purpose it is conveni- 
ent to express I"! in terms of the leading- and trailing-edge singularity strengths. 

Let us introduce the ratio, £, of T'\))(¢) to the other part of the circulation, i.e. 


dé \*dé) 
It is easily shown a 
y(1-,t) _ trailing-edge singularity strength _ (69) 
y(0+,t) leading-edge singularity strength 


so that / is a measure of the relative singularity strengths. £ = 0 corresponds 
to a trailing-edge Kutta condition and # = 1 corresponds to a leading-edge 
Kutta condition. It seems reasonable to restrict # to the range O < # < 1, since 
f outside this range corresponds to ‘supercirculation.’ This will be done here 
although it is not a critical assumption for the analysis. 


Harmonic oscillation 


Consider an oscillating airfoil with the upwash as given by equation (31). Then 
the quasi-steady quantities and wake distributions can again be represented as 
in equations (32) to (35). Using these, the analogous Wagner integral equations, 
(15) and (17), are treated in a manner similar to that leading to (36) and (38). The 
difference is that the integration over €, is forward of the airfoil; this has the net 
effect of changing the sign of the K, terms as indicated in the definitions (37). 
These results can be reduced to 


26GY = 28) KD + K®)] + gof( 28) )K K®), (70) 


— 240 — GP = 4) — 2), (71) 
where £G) corresponds to the quasi-steady contribution of equation (68), 
1 dé 


and the function D (see (40)) is modified to include the wake effects associated 
with /, giving 

D,(w, = [(1—2f) + ay + KY]. (73) 
The lift and pitching moment can be given in a form similar to equations (43) and 


(45). This will not be carried out in detail here; rather, let us turn our attention 
to the stability of the flow. 


Stability of the flow 
For a given airfoil motion and a given /, the wake intensities g, and g, are given 
by the solution of equations (70) and (71), provided the determinant of the co- 
efficients does not vanish. If the determinant vanishes, then there is an eigen- 
solution, i.e. there is a solution with g(x) = 0 but g, . + 0. Let us look for such 
solutions. 
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The determinant of the coefficients of g, and g, in equations (70) and (71) is 


28) Ky) + Ky (1-28) KP +4 exp [tio(1+ 


+(1- =. x I 
1+f)(%+1)]. 
4 
Complex i ln (2n(2n+1)} p=! 
w-plane 
/ 
/ Aiw <0 
/ 
/ 
/ 
i] 
\ 
\ 
Location chiw >0 
of zero 
f=0 
ln (27(2n+1)} —— 


Qnn+a 


Figure 3. Location of an eigen-value, corresponding to a given positive integer », in the 
complex w-plane as a function of « and #. # = 0 corresponds to a trailing-edge Kutta 
condition; # = 1 to a Kutta condition applied at the leading-edge. 


where the asymptotic expansion is valid for |argw| < 7. The zeros of the deter- 
minant for large frequencies can be found by equating the largest terms in the 
expansion. If 0 < £ < 1, then we take 


—B(a+1)+(1—£) (a 


=O 
which gives zeros at 
—1)/P(a+1)}, (75) 


Ro Riw = In{(1—f)(z 


where 7 is a large positive integer. Continuing this analysis, it is found that a 
given pole will move as shown in figure 3, as # goes from 0 to 1. A similar analysis 
shows that there is another set of poles which corresponds to reflexion about the 
imaginary w-axis in figure 3. 

Therefore the determinant of equation (74) has infinitely many zeros, with 
corresponding eigen-solutions. For Ziv > 0, the solutions diverge exponentially. 


Thus, whenever 
A(a+1)/(1—-f)(a—-1) < 1, (76) 
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FIGURE 4(a). The modified Theodorsen function, T(3i, a), for vertical oscillations, case I. 
This is typical of the general case. (b) The modified Theodorsen function, T'(40, a), for the 
case of a downwash distribution linear about the airfoil quarter-chord, case III. This is 
the special case in which the function is not asymptotic to a circle. 
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is Figure 5(a). The function S(jiw,a) for vertical oscillation, case I. (b) The function 
S(}iw, «) for a downwash distribution linear about the airfoil quarter-chord, case ITT. 
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there will be diverging eigen-solutions and the flow will be unstable. Using (69), 
this criterion can be written ast 


Trailing edge singularity strength Rear wake speed | 
Leading edge singularity strength Forward wake speed | 


We conclude, from this brief study of flow with « > 1, that steady flow must 
be unstable when the criterion stated in (76) and (77) is satisfied. It appears that 
under these conditions the vortices that move forward and rearward are capable 
of interacting in such a way as to induce a diverging oscillatory circulation about 
the stationary airfoil. Presumably such oscillation would continue until violent 
separation (stalling) occurs, invalidating the assumptions of the present theory. 

On the other hand, it is conceivable that the resultant effect of viscosity 
and electrical resistance in a real fluid, when a? > 1, is to enforce an effective 
Kutta condition at the leading, rather than at the trailing, edge. This possibility 
is suggested by the studies of Lary (1960), Lewellen (1959), Greenspan & Carrier 
(1959), and Hasimoto (1959) (all of which pertain, however, to steady flow). 
If this occurs, criterion (76) may not be satisfied, and we cannot conclude here 
that the flow is unstable. It seems clear that further study in this area must 
await clarification of the roles of viscous and magnetic boundary layers and sepa- 
ration in fixing the circulation about cylindrical bodies. 


9. Conclusions 

The main conclusion to be drawn from this work is that unsteady effects 
become more important as the magnetic field strength increases. The steady- 
flow theories of Lary (1960) and of Sears & Resler (1959) predict small forces for 
a” near unity; the present theory predicts large forces due to unsteady effects. 
In fact, figures 4 and 5 show that at moderate frequencies the forces may be 
larger than in the corresponding non-magnetic case. 

For a > 1, there is a definite possibility that the flow will be unstable, i.e., 
that no steady solution exists. In particular, if the boundary layers (magnetic 
and/or viscous) act in such a way as to remove the trailing-edge singularity, the 
flow will be unstable. 
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Notes on approximations for turbulent flows with 
large temperature differences 


By D. J. TRITTON 


Department of Aeronautical Engineering, Indian Institute of Science, Bangalore 
(Received 23 August 1960 and in revised form 5 May 1961) 


(i) The diffusion terms in the mean velocity and temperature equations of 
turbulent flow are analysed to decide when variations of fluid properties can 
produce appreciable errors. 

(ii) A theoretical demonstration is given that in the mean-flow continuity 
equation for a gas the error in assuming constant density is small if the flow is 
turbulent, even when the temperature variations are large. 

(iii) Separate discussion is given of the case of local heat sources in turbulence, 
as large errors can occur there. 


1. Introduction 


Intheoretical work on low-speed flows with temperature variations it iscommon 
to suppose that all temperature differences are small enough for changes in fluid 
properties to be neglected. In experimental work the supposition is frequently 
not fulfilled. In gas flows, in particular, temperature differences may well be 
insufficiently small compared with the absolute temperature for density varia- 
tions to have no effect. In any fluid, variations of viscosity and thermal con- 
ductivity may be sizable. These notes consider, for turbulent flows, some aspects 
of the errors so caused. They do not give corrections—the complexity of the 
accurate equations precludes that—but discuss the magnitude of the errors for 
cases in which some analysis is necessary before this can be properly estimated. 
It might be thought, for instance, that for mean flows, where the direct action of 
molecular effects is negligible, temperature variation of viscosity would not 
matter, but this is seen in § 2 to be fallacious. Then, in §$3 and 4, attention is 
drawn to a respect in which theory fails less than would appear at first sight and 
also to a situation in which there is particular danger of large errors. 

The notation is standard. An overbar (~) indicates a mean value and (’) the 
difference between the instantaneous value and the mean value. 

The estimation of errors is just a matter of writing the relative sizes of retained 
and neglected terms in terms of temperature, velocity, and length scales. During 
this procedure, it must be noted that there are two temperature scales (and, 
correspondingly, two density scales); the scale of 7 is that of the absolute tem- 
perature 7), whereas that of 7”, or of 7’ when it appears as a gradient, is A7’, the 
difference between extremes of temperature. Strictly, a further distinction 
should be made between these last two and, similarly, between the scales of mean 
and fluctuating velocities, but this is left out of account; otherwise it is necessary 
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to present the analysis somewhat differently for different flows; the distinctions 
are unimportant in the end. The use of a single velocity scale, U, further does 
not allow for flows such as a wake in which the mean velocity difference has a 
smaller scale than the mean velocity, but it can be shown that the situations in 
which the results are different are of no physical importance. For length scales 
we distinguish between longitudinal gradients having scale L and transverse 
gradients having scale 0, as the results are of interest mainly for boundary-layer- 
type flows. (These appiy, of course, to gradients of mean quantities; we shall be 
meeting a turbulence length scale in § 2.) 


2. Mean velocity and temperature equations with variable fluid properties 


Consider the term 


A = 
Ox? AT 0x; 02; 
B C 


C will behave, at least so far as the subsequent order of magnitude analysis is 
concerned, similarly to 
ere eu; 


—/, where a= 


a function of 7’). 


T du 
pal 
ou; _ OP cu; oT” 


Now 74. — 


and the length scale of the second part is a canake not a mean flow, one. 
Consequently A can be much larger when y is variable than when it is constant. 
In fact C ~ ATU/T,A?, where A is the dissipation-length parameter. [There 
two assumptions involved in this stage. First, no distinction is being made 
between the kinetic energy and 7’? dissipation-length parameters. Hence the 
argument as it stands does not cover the case of Prandtl number very different 
from 1, though it could be modified to do so (nor is the configuration considered 
in § 4 covered by this). Secondly, it is assumed that the correlation 


eT” Cu; eT" 
ox, \\ex,) \ex; 


(summation convention not applying here) is appreciable for at least some 
components. The condition for it to be so may be inferred from the equation 


y 


ct CX; Cx; OX; OX; 


Cx; 


aT vT" ae (+ further terms if buoyancy forces are significant) = 0. 
i 


(1) 
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(That this is the constant fluid property form does not matter for the present 
considerations.) This represents the balance of a component of the temperature 
flux (see Ellison 1957, § 6 and, in particular, equations (16) and (19) for an example 
of this). The seventh and eighth terms indicate molecular effects on this balance; 
they are at least partly dissipative, the second, third, and perhaps the buoyancy 
force terms being the corresponding productive ones. Now these two molecular 
terms are closely related to the correlation under consideration, at least when the 
Prandtl number is near 1, because 


a2 Amy 7? 
oF - oT" Cu; 
tT’ = (uT’)-2— —. 

Oar Cx? Cx? Cx; C2; 


It may be inferred that 2 is appreciable when the balance of w;7" plays a role in 
the turbulence comparable with those of kinetic energy and 72. Thus it is 
appreciable in shear flows with heat transfer but not in grid turbulence with 
uniform mean temperature. | 

Continuing the estimation of the magnitude of C, a little manipulation of the 
various scales in turbulence (Hinze 1959, pp. 185-6), or a consideration of the 


energy balance in a shear flow, gives A? ~ dv/U. Hence 


apy ATU? 


Cx 


We compare this with the magnitude of another term in the same equation, 
ple(u;u') /ex;] (D, say): D ~ py U?/d; hence C/D ~ aAT'/T). Thus for large tempera- 
ture variations in a fluid having « ~ 1, it is incorrect to suppose that molecular 
terms in a mean equation for turbulent flow are small to the order v/v, (where 
Vp is the eddy viscosity), as is the case when fluid properties are constant. 

The reader may well be wondering why this discussion has been given in terms 
of variable viscosity, when it wouid seem that variation of thermal conductivity 


is a simpler case; for 


and the second term is analogous to the one considered above without the com- 
plication of the correlation. However, we shall now see that the same effect does 
not oceur here. We may make the transformation 


a ((Tk 


(kyis the value of k at the reference temperature 7), showing that a term involving 
turbulence length scales need not be brought in and that Z is of the same order of 
magnitude as when k is constant. 

This conclusion is disguised in (2) by the correlation of k-fluctuations with 


(dk/dT’) (0T |cx,)*. This in turn suggests that, if any other temperature-dependent 
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quantity isa coefficient of ¢?7'/¢2?, care is again needed. Consider in this connexion 
the full mean temperature equation 


U;,—+u;— 4 
“Cx, ‘Cx, pC, cx,\ Oz; (4) 


We may anticipate from the above that, although the temperature variafion of 
k does not produce appreciable terms on the right-hand side, that of p does. It is 
seen that this is indeed the case by making the transformation 


d(1/pC,) or __ & 
pC, (ta) Ga) ( = 
Analysis as before shows that the second part of this is small compared with the 
left-hand side of (4) only to the order A7'/T, when £ = —pC, T[{d(1/pC,)/dT}] is 
of order 1 (f = 1 for a perfect gas). 
It is probably misleading to regard this last result as an occasion when molecular 


effects are important in turbulent flow, for the same error in the usual approxima- 
tion may be detected by retaining pC, on the left-hand side of the equation, thus: 


oT 66 (, eT 


There is then an approximation of order Ap/p, in taking the left-hand side of 
(5) as poC,, CT + pyC, ui(OT"| but very little approximation in dropping 
the right-hand side. The important point is just that it is fallacious to argue that 
the equation can be put in the form (4) and the right-hand side then dropped 
because it is a ‘molecular term’. 

We now return to the velocity equation. The above comments about 
p-variability apply again. But we have already seen that ”-variability also has an 
effect; there is no analogue to (3). Intuitively it seems probable that -variability 
will not come in if the boundary conditions on velocity and temperature are 
similar and if, further, k(7’) and ”(7’) have similar functional forms. In general, 
however, it does come in. It is particularly important to remember this for fluids 
for which the variation of ” is more marked than that of p (« large compared 
with /); water is an important example in which this is strongly so. Then there 
may be an approximation in the velocity equation but little in the temperature 
equation. 

All the preceding discussion has been concerned with the approximation in 
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3. Continuity equation of mean flow 


The continuity equation of mean flow is usually taken in the form 


It will now be shown that this is a good approximation in turbulent flow of gases 
even when the temperature variations are large—far better than is immediately 
apparent. Physically this means that expansion effects shift the mean flow 
streamlines less than in laminar flow (for which ou,/ox; = 0 is quite a poor 
approximation). The result is probably most useful for the case of a boundary 
layer on a heated wall—the temperature difference can remain large at all 
distances downstream—but the treatment here will be general. 
The accurate continuity equation may be written 
ou, 1 Dp 


0x, p Dt (6) 


The second term is similar to the left-hand side of the temperature equation 


and may thus be replaced by one similar to the right-hand side. This is a molecular 
conductivity term and the result of this section derives in essence from the 
unimportance of molecular effects in turbulent flow (it will be seen that the 
exception to this described in $2 does not usually arise here). Substitution of 


pT = poT') m (7) gives 
1Dp___ 
pm C, Oat, 
and so from (6) 
1 0 oT (ft 
Ox; Cy Ox; (« (| kar) 


We then obtain, on taking average values, 
OU 2 
G 
We now consider the order of magnitude of either side of (9): 
FrU/L; G~ «xAT/Te. 
(It is possible to think of a form of k(7') for which the use of AT’ as the scale of the 


difference of | = dT is unsatisfactory, but such is most unlikely to apply to areal 


fluid. This is, for instance, the correct scale for any power law.) Thus we have 
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G/F ~ ATKL/T, Ue. But L*/62 ~ UL/vp. (It is because this increase in the 
boundary-layer thickness through an eddy viscosity replacing a molecular one is 
essential to the argument that it does not also apply to laminar flow.) Therefore, 


G AT 

Pro,’ 
which will be small even if A7'/J, is not. The fractional error in putting F = 0, i.e. 
0u,;/0x,; = 0, will be correspondingly small. 

The preceding analysis made use of the perfect-gas equation. This should apply 
in most cases of practical importance. However, brief mention may be made of 
what happens if a different equation of state applies; the situation is somewhat 
more complicated, and it may clarify the nature of the analysis to see why. For 
any other equation of state, the equation corresponding to (8) has temperature- 
dependent quantities in the coefficient on the right-hand side. Consequently 
averaging cannot be carried through the differentiation and, in the same way as 
in §2, a term involving turbulent length scales comes in. The magnitude of this 
may be estimated as before;+ it is of order (A7'/7))? smaller than F’, which still 
compares favourably with a fractional error of A7'/T, in laminar flow. The still 
smaller error in the perfect-gas case results from a coincidental cancelling of two 
non-linearities, that in the equation of state and that due to the temperature 
dependence of p in the term pC,,(DT'/ Dt) of the temperature equation. 

Expansion effects may remain important in a laminar sublayer; indeed the 
above argument suggests that, in view of the large gradients, they may be 
particularly serious there even for comparatively small A7'/7,. However, 
qualitative prediction of their behaviour might be possible by analogy with flows 
with injection through the wall. 


4. Temperature difference introduced locally into pre-existing turbu- 
lence 

So far it has been supposed that the length scales of the mean velocity and 
temperature fields are the same. This is not so when, as in turbulent-diffusion 
experiments, the heat source is a small element placed in an already turbulent 
flow field. A full analysis for this case of all the aspects considered above is com- 
plicated and, so far as I can see, not very rewarding. However, this is a situation 
in which particularly large errors can arise. We therefore consider one example 
of the way this can happen, by examining the relative sizes of u;(0u;/0x,) (H, say) 
and 6(u;w;)/cx,; (J, say), the former being the term neglected in the usual practice 
of taking the latter as the gradient of the stress. If the one is not small compared 
with the other, the supposition that the velocity field is unchanged by the intro- 
duction of the heat source is without justification (although the ratio H/J may 
overestimate the change; as the temperature difference is increased, the flow 
might adjust itself so that the gradient of the true stress remains constant rather 
than the gradient of w;u}). 


+ Though this is a point at which the distinction between the scales of mean and 
fluctuating quantities (see §1) 7s necessary to the argument. 
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This particular example was chosen primarily because it is probably the most 
serious one, but also because it illustrates the way in which the unapproximated 
continuity equation (6), with the division into mean and fluctuating parts applied 
to both velocity and density, can be used for estimating the magnitude of such 
effects. Equation (6) may be written 


where r = log p/p) (and thus has Ap/p, as its scale and a spatial distribution 


similar to that of temperature). Subtraction of the average of this equation from 
the unaveraged form gives 


a * “On, ‘On 

Ou; _ ,0 or 

and so Ue = — UW — — — — — —. 10 

H K L M N 


This can be used for estimating the size of H. The difference from previous cases is 
brought about by the heat wake having different length scales from the main 
flow; we will call its width 6,. From (10) we get H ~ ApU?/p,6,,. This is im- 
mediately indicated by the term Z.} Furthermore, J ~ U?/d, and so 


which is not necessarily small. As an illustration of the possible effect in an 
experimental set-up, consider 4/0, in the region of measurement in Corrsin & 
Uberoi’s (1951) work on local heat sources. It is about 10. These authors are not 
explicit about the magnitude of Ap/p,, but it seems likely from the temperature 
and size of the source that it was somewhere between 0-02 and 0-1. Thus we see 
that H/J can be quite appreciable. 

A similar analysis may be carried through for the effect of non-zero cw;/dx,; on 
the use of 0(w;7")/0x, in the temperature equation. However, this is as usual small 
to the order Ap/py; 6;, is the length scale in all the terms. Diffusion of heat from the 
source is seriously affected only through changes in the mean velocity field. 


I thank Dr A. A. Townsend and a referee for comments that have greatly 
helped with the development of these ideas. My work at Bangalore is made 
possible by the award by the Royal Society, to which I am most grateful, of a 
Autherford Scholarship. 
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y The other terms are not so readily estimated, but there is no reason to expect either 
that they are any larger or that they just cancel Z so as to make H smaller. In point of 
fact, arguments invoking the similarity of these terms to ones appearing in equation (1) 
suggest that K, and N are likely to be of the same order as L. 
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The transverse force on a spinning sphere moving 
in a viscous fluid 


By S. I. RUBINOW 
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The flow about a spinning sphere moving in a viscous fluid is calculated for small 
values of the Reynolds number. With this solution the force and torque on the 
sphere are computed. It is found that in addition to the drag force determined 
by Stokes, the sphere experiences a force F;, orthogonal to its direction of motion. 
This force is given by F,, = x V[1+O(R)]. 
Here a is the radius of the sphere, Q is its angular velocity, V is its velocity, p is 
the fluid density and F& is the Reynolds number, R = px-'Va. For smail values 
of R, the transverse force is independent of the viscosity . This force is in such 
a direction as to account for the curving of a pitched baseball, the long range of 
a spinning golf ball, etc. It is used as a basis for the discussion of the flow of a 
suspension of spheres through a tube. 

The calculation involves the Stokes and Oseen expansions. A representation of 
solutions of the Oseen equations in terms of two scalar functions is also presented. 


1. Introduction 

We shall show that a spinning sphere moving in a viscous fluid experiences a 
force orthogonal to its direction of motion, which we call a lift force. It is also 
orthogonal to the spin axis and therefore accounts for the curving of a pitched 
baseball, the long range of a spinning golf ball, etc. Since we shall calculate the 
lift for small values of the Reynolds number, our result will not apply to these 
phenomena, which occur for large Reynolds numbers. Our result applies either 
to the slow motion of a small sphere in a fluid, such as a solid particle in a suspen- 
sion, or to the motion of a satellite or other astronomical object in a gas of low 
density. We shall also show that the spin does not effect the drag force and that 
there is no correction to the retarding torque to the order we consider. Using 
these results we shall determine the motion of a sphere with given initial linear 
and angular velocities. We shall also consider the flow of a suspension of spheres 
through a tube, in the light of our results. 

Our results are obtained by solving the Navier-Stokes equations for the motion 
of the fluid around the sphere. We determine the first few terms, in the expansion 
in terms of the Reynolds number, of the solution. This expansion consists of two 
parts, the Stokes and Oseen expansions, which were introduced and used by 
Lagerstrom & Cole (1955), and Proudman & Pearson (1957). Our method of 
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solution is similar to theirs with the main difference being that we obtain directly 
the expansion for the velocity and the pressure rather than for the stream func- 
tion. This is necessary because our flow is asymmetric and therefore there is no 
stream function which describes it. Those authors, as well as Stokes (1851) and 
Oseen (1910, 1913) considered the motion of a non-spinning sphere, about which 
the flow is axially symmetric and describable by a stream function. The present 
example appears to be the first asymmetric three-dimensional case for which the 
flow at small Reynolds numbers has been found by using these expansions. In 
the appendix we prove that it is possible to represent any solution of the Oseen 
equations in terms of two scalar functions. 


2. Formulation 


We wish to consider the motion of a sphere of radius a, velocity V and angular 
velocity Q through an incompressible viscous fluid of density p and viscosity 
coefficient “. The fluid is assumed to be at rest far from the sphere. It is con- 
venient to introduce a co-ordinate system with its origin at the centre of the 
sphere, with the negative x’-axis pointing in the direction of V and with the 
(x’, y’)-plane containing Q. In this co-ordinate system the flow is steady and the 
fluid has the velocity (V’, 0, 0) at infinity while the velocity of the surface of the 
sphere 7’ = ais Q xr’. We denote the velocity in the fluid by uw’ and the pressure 
in the fluid by p’, taking the pressure at infinity to be zero. Now u’ and p’ satisfy 
the time independent Navier-Stokes equations, the equation of continuity, the 
condition of no slip at the sphere and the appropriate conditions at infinity. 
These equations and conditions are 


pAu’—Vp’' = p(u’.V)u, (1) 

= (2) 

at r <a, (3) 

u’'= (V,0,0) at =o, (4) 

and at (5) 


Let us introduce the dimensionless Reynolds number RF as well as the dimen- 
sionless quantities u, p. w, x, y, z and r, defined by 


u=V4u’, w=aVQ, R= puVa.) 6) 
Upon introducing these variables into (1)-(5), we obtain 
Au-—Vp = R(u.V)u, (7) 
V.u=0, (8) 
u=oxr at r=1, (9) 
u=(1,0,0) at r=0, (10) 
and p=0 at r=0. (11) 


We seek the solution u and p of (7)-(11) as expansions in RF valid for small values 
of R. 
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3. Stokes expansion 
The Stokes expansion of the solution is of the form 
P = Pot Rp, +o(R). . (13) 
Further terms beyond those shown involve powers of log R in addition to powers 
of R, as has been demonstrated by Proudman & Pearson (1957) for the non- 
rotating sphere. However, only the terms shown in (12) and (13) will be needed 


here. We now insert (12) and (13) into (7)-(9) and equate coefficients of 2° in 
each equation, obtaining 


Auy— Vp = 9, G4) 
V.u, = 0, (15) 
and at (16) 


From the coefficients of R!, we obtain 


Au, — Vp, = (Up. V) Up, (17) 
V.u, = 0, (18) 
and u=0 at (19) 


We have not required of the Stokes expansion that it satisfy (10) and (11), 
the conditions at infinity, because this expansion is not uniformly valid in the 
neighbourhood of infinity. Therefore conditions at infinity on the individual 
terms cannot be obtained in the same way. Instead we must obtain another 
expansion, the Oseen expansion, which is valid at infinity. Then by matching 
the two expansions the necessary conditions will be obtained. To obtain the 
Oseen expansion we must first introduce stretched variables. 


4. Stretched variables 
Let us introduce the new stretched variables X, Y, Z and s and the new 
functions U and P, defined by 
A=Re, Y=Ry, Z=Rz, = Rr, 
U(X, Y,Z, R) = u(R-1X, R-'Y, R-1Z, R), (20) 
P(X, Y,Z, R) = R-'p(R-1X, RY, R-Z, R). 


In terms of these variables, equations (7)-(11) become 


AU-VP = (U.V)U, (21) 

V.U =0, (22) 

U=Rwxs at s=R, (23) 

U=(1,0,0) at (24) 

and P=0 at 00. (25) 


29 Fluid Mech. 11 


y 
lar 
ity 
on- 
the 
the 
the 
ure 
sfy 
che 
ty. 
(1) 
(2) 
« 
3 
(3) = 
(4) 
6) 
- 
8) 
9) 
0 
i! ) 


450 S. I. Rubinow and Joseph B. Keller 


5. Oseen expansion 
The Oseen expansion of the solution U and P of (21)-(25) is of the form 
P= (27) 


Upon inserting (26) and (27) into (21), (22), (24) and (25) and equating coefficients 
of R®, we obtain 


AU, — = (Uy. V) (28) 

V.U, = 0, (29) 

U, = (1,0,0) at s=00, (30) 

and at 8 = 00. (31) 


Kquating coefficients of R! yields the equations 


AU, —VP, = (Uy. V) U, + (U,.V) Ug, (32) 

V.U, =0, (33) 

U, = (0,0,0) at (34) 

and P,=0 at 8s = 00. (35) 


We have not inserted the Oseen expansion into the condition (23) on the sphere 
because we do not expect this expansion to be uniformly valid there. Con- 
sequently this condition must be replaced by another condition. This latter 
condition is that the Stokes and Oseen expansions must match since they are 
both expansions of the same solution expressed in different variables. This 
matching will be explained later. 

We note that equations (28)-(31) satisfied by U, and P, are the same as the 
corresponding equations for U and P. Therefore we cannot expect to solve them 
more readily than we could solve the original equations. However, we observe 
that a particular solution for Uy and P, is 


U, = (1,0, 0), (36) 
P,=0. (37) 


The Oseen expansion is based upon choosing this particular solution for U, and 
P,. With this choice (32) simplifies to 
AU, -VP, = 


ax 


Before determining U, and P,, let us explain the principle of matching. 


6. Matching 


The matching principle is that two different asymptotic expansions of a given 
function must be asymptotically equal in their common domain of validity, if 
any. We assume that the Stokes expansion is valid from the sphere out to some 
large distance. We also assume that the Oseen expansion is valid from infinity 
in to some small radius in the stretched variables. However, this radius is 
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considered to be large in the unstretched variables. Then the common domain 
of validity is a spherical shell within which the unstretched radius is large but 
the stretched radius is small. Therefore if the Stokes expansion is expanded for 
large values of the radius r while the Oseen expansion is expanded for small 
values of the stretched radius s, the resulting expansions must be asymptotically 
equal. This principle is utilized by alternately determining terms in the two 
expansions and matching them to the previous terms. 

As a first result of applying this principle we can determine the behaviour of 
the zeroth-order Stokes approximation, uy, and po, for large values of 7. Since 
U, and # are constant, the result is easily seen to be 


Uy = Upto(l) as r>o, (39) 


and Py = 0(1) as roo. (40) 


7. Zeroth-order Stokes approximation 

The solution u, and p, of (14)—(16), (39) and (40) is unique and can be found by 
the method explained by Lamb (1945, pp. 595, 596). Since the equations are all 
linear, the solution is just the sum of the known solution for a uniform flow past 
a non-rotating sphere and the known flow produced by a rotating sphere in a 
fluid at rest at infinity. Thus we have 


3 1 3zf/1 1 1 
3a 
Po = — (42) 


8. First-order Oseen approximation 

To determine U, and P, we must solve (38), (33)—(35) and employ the matching 
condition. We begin by taking the divergence of (38) and using (33), from which 
it follows that AP, = 0. Then, following Lamb (1945, p. 610), we introduce a scalar 
@ and a vector W and write 


% 
U, = —Vé+W. (44) 


Since P, is harmonic we require that ¢ also be harmonic, i.e. 
Ag = 0. (45) 


Then (38) and (33) become two equations for W: 


V.W=0. (47) 


In the appendix we prove the following theorem. 
29-2 
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Theorem. If W satisfies (46) and (47) then there exist two scalar functions 


y and y such that = (1,0, 0) x Vyr+V x [(1, 0, 0) x Vx] (48) 
0 
and Ay- OX = 0, (49) 
Ay- 0 (50) 
Let us define y’, which also satisfies (50), by 
xX (51) 
Then (48) can be written as ox 
W=VvVx +(-x.- 35.54): (52) 


Those special solutions of the form (48) or (52) for which y = 0 were given by 
Lamb (1945, p. 611). 

To construct U, and P, let us make the following choices of ¢, y and y’, which 
must be solutions of (45), (49) and (50) respectively: 


Als, (53) 
y=9, (54) 
and x’ = (B/s)exp[4(X —s)]. (55) 


We now compute U, and P, from (43) and (44), using (52)-(55) and employ the 
matching condition. This yields A = B = 3. Thus U, and P, are found to be 


3 1 2 


(56) 
and (57) 


This solution satisfies the condition at infinity (34) and (35). If any other choice 
were made for ¢, iv and y’, the matching condition or the conditions at infinity 
would have been violated. 


9. First-order Stokes approximation 

To determine u, and p, we must solve (17)-(19) and employ the matching 
principle. We shall first find a particular solution U4) and p49 of (17) and (18). 
Then we shall add a solution u,, and p,, of (18) and the homogeneous form of (17) 
in order to satisfy the boundary condition (19). Finally we shall show that the 
sum of the two solutions so obtained satisfies the matching condition so no other 
solution need be added. To begin, we compute the inhomogeneous term (uy. V) Uy 
which occurs in (17). By using (41) for uy we obtain 


16 
12 12 20 24 4 


4 


1 
(0,0, 1) + x (w xr). 


(58) 


Here f is the angle between 2 and Up). 
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Let us now take the divergence of (17) and make use of (18) and (58). Then 
we find that 
9 ‘46 2 
Ap, = —V.[(u,.V) ug] = | 


r 


( . (59) 


A particular solution of (59) is 


_ 3 6 2 1 3 wzsnPp (3 2 (w.r)? 


(60) 
We now insert (60) into (17) and obtain for u, the equation 
Au, = Vpi9 + (Ug. V) Uo. (61) 
A particular solution of (61) is 
3 4 3 1 & 
wsnf[S 1 wz sin 2 
— 
16 E pt 16 E 4 
wx(wxr) w* (w.r)? 


By direct computation we find that uj, satisfies (18). 
Next we find u,, and p,, satisfying the homogeneous form of (17) and (18) and 
such that U4 )+U,, satisfies (19). Thus we require that 


By the method of Lamb (1945, pp. 595, 596), we find that 


x sin wzsin 7 


@x(wxr) 3 (@w.r)?—3 5 
47° 8 (3+ 8 yt 


(; *) wzsinp ow 3(w.9r)? 


Upon adding uj, to u,, and pj, to p,,, we obtain u, and p,. (There is no need 
here to write out the resulting expressions for u, and p,, but we shall presently 
refer to them as (66) and (67), respectively.) By employing the matching pro- 
cedure we find that this solution matches the Oseen solution, and therefore it 


is not necessary to add anything else to it. 
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10. Lift, drag and torque 


We shall now use the two terms of the Stokes expansion which we have found 
to compute the force and torque on the sphere. To do so we must first determine 
f, the force per unit area exerted by the fluid on the surface of the sphere. Lamb 
(1945, p. 596) gives an expression for f which becomes in dimensionless variables 


Vel | 1 
Upon using (41), (42), (66) and (67) in (68) and integrating (68) over the surface 
of the sphere we obtain the total force F which is given in terms of the velocity V 
of the sphere by 


F = x V+o(anV R). (69) 


The component of F in the direction — V is the drag force F,, and the component 
normal to V is the lift foree F;. From (69) and the definition (6) of R, these 


and F, = 7a°pQ x V[1+O(R)]. (71) 


It is interesting to note that the leading term in F,, is independent of and is 
similar to the lift formula for two-dimensional potential flow about an airfoil. 
By integrating r x f over the surface of the sphere we obtain for the torque T 


on the sphere the result T = —8npo°Q[1 +0(R)]. (72) 


The term — 6z7asV in F, was calculated by Stokes (1861) and the second term 
—{mpaVR by Oseen (1910, 1913). The expression for F,, which is new, is the 
main result of this paper. The term — 87ya°Q in T was obtained by Kirchhoff 
(1876). Our results show that the spin causes no correction to the drag of order 
R and that there is no correction to the torque of order R. 


11. The motion of a spinning sphere 


Let us now consider the motion of a sphere with initial velocity V, and initial 
angular velocity 92) in a viscous fluid initially at rest. Let 17 denote the mass of 
the sphere, J its moment of inertia about any diameter, V(t) its linear velocity 
and 92(¢) its angular velocity about a diameter at time ¢. We assume that the 
force and torque on the sphere are given by (69) and (72), ignoring the terms 
o(ayV R) in (69) and o(R) in (72). Then the equations of motion of the sphere are 


M = — V (73) 
and I = — (74) 


The solution of (74) is 


Q(t) = Qyexp (75) 
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If V(t) denotes the component of V(¢) parallel to Q), (73) yields 


V(t) = VY exp [ — (76) 
Let V(t) denote the component of V(t) perpendicular to Q). Then (73) yields 
[V2(t)| = |VP| exp [— (1 + (77) 
Now we define the unit vector v(t) by 
= |V(t)| v(t). (78) 
From (73) it follows that v(t) satisfies the equation 
exp Q, x Vv. (79) 
The solution of (79) is, in rectangular components, 
v(t) = [cos A(t), sin A(é)]. (80) 
Here the angle 6(¢) is given by 
O(t) = (: —exp |- tl) (81) 


From these results the path of the sphere can be obtained easily by integration. 


12. Flow of a suspension through a tube 


When a viscous fluid containing suspended particles flows through a circular 
tube, the particles are not uniformly distributed over the cross-section of the 
tube. Observations on blood flow in narrow capillaries indicate that the red 
blood corpuscles concentrate near the axis of the capillary. This concentration 
effect was explained hydrodynamically by supposing that the iluid exerted an 
axially directed transverse force on the particles. However, when the flow was 
determined by solving the linearized Stokes equations, no transverse force was 
found on a sphere in a shear flow, such as the Poiseuille flow in a circular tube.7 
A similar calculation of the force on an ellipsoid in a shear flow by Jeffery (1922) 
also yielded no transverse force. 

To account for the effect, Jeffery proposed a new principle—that a particle 
moves at that distance from the axis which minimizes the rate of energy dissipa- 
tion in the flow. This principle was subsequently rediscovered and generalized 
and has become a cornerstone of irreversible thermodynamics. If correct, 
Jeffery’s principle should be deducible from the Navier-Stokes equations, just 
as all thermodynamic principles are deducible from more detailed theories. To 
apply it to the present problem, Jeffery utilized the result of Einstein (1906) for 
the excess dissipation produced by the introduction of a sphere into a uniform 
shear flow. Einstein has shown that the excess dissipation is proportional to 

+ This result can be obtained without calculation by decomposing the incident shear 
flow into two parts, one symmetric and the other antisymmetric about any plane through 
the centre of the sphere and containing the direction of the incident flow through the 
centre. By symmetry, neither of these flows can yield a force on the sphere normal to 
the plane. 
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the square of the vorticity of the incident flow at the particle. Since the velocity 
distribution in the Poiseuille flow is parabolic, the vorticity is proportional to 
the distance from the axis and the excess energy dissipation to the square of the 
distance from the axis. Thus Jeffery concluded that a particle would move along 
the axis. 

This same conclusion was arrived at by Tollert (1954) and Saffman (1956). 
Tollert observed that, according to Einstein, a sphere in a shear flow u rotates 
with the angular velocity 2 = IV x u. Furthermore, according to the calculation 
of Simha (1936), a sphere of radius a in a tube of radius R, lags behind the flow 
with the relative velocity V = — §(a?/R§) Up, where U, is the flow velocity at the 
centre of the tube. Then Tollert postulated a transverse force of the form (71) 
with another factor instead of 37. In a Poiseuille flow, V x u = —2b x U,/R?, 
where b denotes the radial vector from the axis. If the correct formula (71) for 
the transverse force on a sphere at radial position b were used, Tollert’s argument 


would yield the result F,, = —2a°pU3b/3Ri. (82) 


Saffman derived a similar result by approximate solution of the flow equations, 
obtaining the numerical factor 7-7(37) instead of 37. The result (82) yields a 
force toward the axis which vanishes on the axis, in agreement with Jeffery’s 
conclusion. 

Recently Segre & Silberberg (1961) measured very accurately the distribution 
of spherical particles in a fluid flowing in a circular tube. They found that the 
particles do not concentrate along the axis. Instead they concentrate at the 
radius 6) equal to about $y. Thus the previous observations were inaccurate, 
and the theories of Jeffery, Tollert and Saffman are incomplete. Consequently 
the formula (82) is incorrect, or at best incomplete. This formula and the theories 
referred to only account for the inward motion of spherical particles lying outside 
the radius by. 

In order to understand why (82) fails, we must recognize that it is based upon 
(71), which applies to a sphere in a uniform flow. But the Poiseuille flow in a 
tube is not uniform. It is a shear flow with a parabolic velocity profile. The 
gradient of this profile accounts for the spin of the sphere and the curvature 
of the profile accounts for the lag. However, these are evidently not the only 
ways in which the non-uniformity of the flow affects the sphere. We therefore 
attempted to calculate the flow about a spinning sphere in a shear flow with 
a parabolic velocity profile, and from it to compute the force and torque 
on the sphere. We followed the procedure employed in the preceding sections 
and had no difficulty in obtaining the zero-order Stokes solution. We began to 
calculate the first-order Stokes solution and obtained some terms in it, but we 
did not complete the calculation because the labour became prohibitive. How- 
ever, from the terms we calculated, we computed the transferse force, which 
came out to be given by (82). But there are many more terms in the force which 
we did not calculate. We expect that these further terms would modify (82) in 
such a way as to account for the observed effect, but this is merely a conjecture. 
This conjecture is strengthened by the fact that (82) yields a force of the correct 
order of magnitude to account for the observed effect, as we shall show below. 
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We have also examined the minimum-dissipation-rate analysis of Jeffery 
and found that it utilizes the linearized flow about a sphere in an unbounded 
uniform shear flow. But to zero order in the Reynolds number, the correct flow 
is the linearized flow about a sphere in a Poiseuille flow in a tube. If this zero- 
order flow were used in calculating the rate of energy dissipation, it might yield 
the correct equilibrium distance of the sphere from the axis. : 

In view of the experimental results, let us now modify (82) so that it gives an 
outward force forb < by. The simplest way to do this is to multiply the right-hand 
side of (82) by (b —b9)/by. With this factor the force is inward for b > by and out- 
ward for b < by. Of course this expression for the force is not correct in detail, but 
it is qualitatively correct. 

We shall now calculate the trajectory of a sphere subject to this force. Its 
radial velocity db/dt can be determined by equating the drag force 62aydb/dt 


to F,. Thus db b(b—by) (83) 


dt 9Rin by 
Now the axial velocity of the sphere dx/dt is practically equal to the flow velocity 
— b?Ro?). Thus we have, upon dividing, 
db atpU,b(b — by) 
We now integrate (84) from 2 = 0, denoting the initial radius of the particle at 
xv = 0 by b,. Thus we obtain 


For 6 nearly equal to by, (85) becomes 
— b?) b, —bo 
8 
Upon solving (86) for 6 and introducing the length x) we obtain 
b—by = (b, — by) exp | — 2/29]. (87) 


Here x,y is defined by 
(1 — Ro*) 


= 


The concentration of particles c(b, 2) in the steady state is easily found, from 
the conservation equation, to be given by 
b,(db,/dt) e(b,, 0) 
b(db/dt) 
By using (83) in (89), setting ¢(),,0) = ¢, and then using (87) to eliminate b, we 
obtain 


c(b,x) = (89) 


bi(d, by) 


{bo t+ (b— by) exp 


c(b, x) 


Cy 


(b)(1—exp [—a/x9]) < b < by 
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Outside the range determined by the extreme initial radii b, = 0 and b, = R,, 
which is given by the inequalities in (90), the concentration is zero. The con- 
centration measured by Segre & Silberberg is in qualitative agreement with 
(90). Furthermore, they have found that their concentration data depend on x 
only through the ratio z/x). where 2» is given by (88). This agreement shows that 
the force is of the order of magnitude given by (82). Since this force is of the first 


order in the Reynolds number, the actual force is also of this order and pre-. 


sumably would be obtained by completing the calculation of the first-order 
Stokes approximation to which we referred above. 


The research in this paper was partially supported by the Office of Naval 
Research. 


Appendix: on the solution of the Oseen equations 


We wish to prove that every solution W of the reduced Oseen equations (46) 
and (47) can be represented in the form (48) in terms of two scalar functions yy 
and y satisfying (49) and (50). To do this we first rewrite (48) in component form 
after setting W = (IV, V3, W3) and obtain 


We shall show that we can determine y and y to satisfy (A 1). First we note that 


if x satisfies (50) then y,,+X.2= Xe—Nzz- Then the first component of (A 1) 
becomes an ordinary differential equation for y which has the solution 


| “Wily. y.2) dy dé +a(y, 2) + 2). (A 2) 
0 0 


Here a and b are two ‘constants’ of integration. 
To determine a and b we apply A, = 0?/éy? + ¢c?/¢ez? to (A2). From (A1) we 
see that A, y = W, while (46) shows that A, W, = W,,—j,,,. Thus (A 2) yields 


x 
W, =e | | — dé + + Agb. (A3) 
0 0 


By evaluating the integral in (A 3) we obtain 

W, = z) —W4,(0, y, z)] + b(y, z) + y. 2) —Wi(0,y, 2). (A 4) 
In order that (A 4) hold for all values of x, y, and z it is necessary that a and b 
Ayaly,2) = W,(0,y,2) (A5) 
and A, b(y, z) = W,(0, y, z) —W,(0, y, 2). (A 6) 
These equations have solutions which are determined up to additive harmonic 
functions. Thus (A 2) determines a class of functions y which satisfy the first 
equation of (A 1). 

Let us now use the second and third components of (A1) to determine y. 


These components yield = —W,-y,, and = W3+y,,.. The integrability 
condition y,,—,, = 0 is then 


W3. + + Woy + Xeyy 0. (A 7) 
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Since y,,,+X.. = W by (A1), (A7) becomes W,, + Wy, + W3, = 0 and this equation 
is satisfied according to (47). Thus we may compute y from the relation 


y 
0 0 
= 0,0)+ | ,0) +X 9, dy 
0 


We have now determined y and yf that satisfy (A 1). We have also verified that 
satisfies (50) since we have shown that = W, and = so 
Ay— x, = 0. Therefore we need only verify that 7 satisfies (49). To do so we first 
observe that y,, = and = Thus = We, — 
Next, from (A 8) we have 


Vy ~ Vex = (2, 0, 0) 0, 0) +’ Worx dy 


—| Ware + (49) 
From (46), W3,.—W5,., = A, W; while from (50) and (A.1) = = Wi, 
Similarly, = AoW, and = Furthermore, from (47) 
Wire = — (Wy, + W5.), and = — (Wh, + W5.),. Then (A 9) becomes 


a 


2 
he— = Vel, 0, 0, 0)+ | — | 


= 0,0) — 0, 0) — 0, 0) + 0, 0) 
+Wy,(x, y,2)— W(x, y,2). (A10) 


Let us now require that y(a, 0,0), which is so far arbitrary, satisfy the equation 
Y,(x, 0, 0) — 0,0) = W3,(x, 0, 0) — 0, 0). (All) 


Then (A 10) shows that y,—y,,. = W3,—W, which we saw above is equal to 
Vuyt Wz: Thus y satisfies (49) and the theorem is proved. 
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Photographic evidence of the formation and growth of 
vorticity behind plates accelerated from rest in still air+ 


By D. PIERCE 


Aerodynamics Department, Royal Aircraft Establishment, Farnborough, 
Hampshire, England 


(Received 5 May 1961) 


The shedding of vorticity from the edges of a plate, which is accelerated normal to 
itself from rest in still air, is investigated experimentally by means of a new flow- 
visualization technique using spark-lighted shadowgraphs and _ heated air, 
benzene or other vapours to introduce sufficient changes in density along the 
surface of discontinuity. 

The photographs show small-scale undulations and the formation of a number 
of centres of vorticity along the well known big-scale vortex sheets which roll up 
along their edges. 


1. Introduction 

Prandtl, in his well known paper of 1904, not only introduced boundary-layer 
theory but devoted the larger part of the discussion to the mechanism and 
consequences of boundary-layer separation. Using the example of the initial 
motion of a fluid as set up by a plate moved from rest at right angles to itself, he 
discussed the development of the shear layer which originates from the plate’s 
edge, and which may be regarded as a Helmholtz surface of discontinuity, or 
vortex sheet, in the case of small viscosity. He pointed out that the vortex sheet 
tends to roll up along its free edge and, further, that such a sheet may be unstable 
in that small disturbances of its shape are likely to be amplified and lead to the 
formation of a series of centres of vorticity along the sheet. In this way, a small- 
scale undulation and vortex formation is superimposed upon a big-scale vortex 
sheet which itself may incorporate a rolled up edge with the formation of a vortex 
core there. 

Prandtl described the flow mechanism and the resulting features from observa- 
tions made in a water channel, and similar observations in water have since been 
made (Wedemeyer 1956). However, most of the later studies have been con- 
cerned with problems of the big-scale flow: the time-dependent development of 
this flow (Anton 1939; Wedemeyer 1956); the corresponding three-dimensional 
flow for the case of the flat triangular wing (Mangler & Smith 1959); and the 
features of the main core (Hall 1961). The small-scale undulations have been most 
extensively studied in the related case of the shear layer round a free jet (Wille 
1960), and it has been noted that the formation of centres of vorticity and ‘spots’ 
in laminar layers may constitute an important phase in the mechanism of 
transition to turbulence. 


+ Crown copyright reserved. 
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The present study returns to Prandtl’s original problem of a plate moved 
normal to itself from rest in still air. It is demonstrated how spark equipmen‘, 
designed for the examination of unsteady supersonic flows, may be used to stud:; 
time-dependent air flows at low speeds by the introduction of small quantities of 
a foreign vapour. 


2. Equipment 

The equipment used can be described in two parts: (a) that which concerns the 
acceleration of the plate from rest, and (b) that which enables the flow to be 
photographed. 

2.1. Plate acceleration 

The mechanism for accelerating the plate is shown in figure 1. The equipment 
consists essentially of three moving parts, with provision for the sudden applica- 
tion of pressure to two of the moving parts by the bursting of a diaphragm. 

The plate under test is screwed to the web of a free piston which initially, by 
gravity, is in contact with a limited-stroke piston, as shown in figure 1. Pressure 
is increased in the compressed air chamber until the polyester film diaphragm 


Retardation cylinder 
L originally filled with 
compressed air, later 
with sponge rubber 


_ Retardation piston 


Free piston carrying 
plate section on 
support web 


N ===) 
/ 
Plate section 
Limited stroke 

piston 

ring 


—=— Crystal for 
SH] triggering spark 
light source 


Bursting diaphragm 


| 

| 

| 

Compressed air | 


FicureE 1. Equipment for accelerating the plate sections. 
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bursts and the pressure is suddenly applied to the limited-stroke piston. This 
piston and the free piston carrying the plate are accelerated together until the 
limited-stroke piston reaches the extent of its travel. The free piston carries the 
plate and continues to move with roughly constant velocity (it is decelerated 
slightly by gravity) until after moving for a few inches it is brought to rest by the 
retardation piston. 

Originally the retardation cylinder was filled with air from the compressed air 
chamber so ensuring identical pressures for acceleration and retardation. Experi- 
ments, however, showed the retardation to be too severe and the air in the 
retardation chamber was replaced by sponge rubber. An O-ring on the limited- 
travel piston acts as a buffer for its retardation and as a seal to prevent the high- 
pressure air from escaping. This method of suddenly applying load to the piston 
obviates the necessity for any breaking link or quick-release mechanism, and it 
has been found that 0-00025in. thick diaphragms consistently give plate 
velocities of 12-13 ft./sec and 0-0005 in. thick diaphragms give 23-24 ft./sec, the 
variation in velocity depending mainly on slight differences in weight of the 
various plate sections. 

The total weight of the moving parts has been kept small and is approximately 
20z. A piezo-electric crystal is mounted in the wall of the chamber between the 
diaphragm and the face of the limited-stroke piston to provide a means of 
triggering the photographic equipment. This acceleration system is clearly 
capable of further development, such as to enable the plate to be accelerated to 
the required velocity in a shorter distance and also to provide for an increased 
length of travel. 

2.2. Flow visualization technique 


The present experiments were undertaken in air at room temperature and 
pressure, i.e. roughly 65 °F and 29in. of mercury, and a shadowgraph system was 
used together with spark photography to show the details of the vortex sheet. 

A spark apparatus which produces an intense light fiash with a duration of 

1 usec was available. This was used as the light source of a shadowgraph system in 
which the parallel light beam was directed over the plate section. The spark flash 
was triggered at an appropriate time by the signal from the piezo-electric crystal, 
the signal being first passed through a pre-set delay unit. [ford H.P.S. photo- 
graphic plates were used and positioned close to the moving plate to reduce any 
effects due to stray light. An ordinary shadowgraph technique cannot be used 
directly because the density changes in the flow are too slight. Thus density 
changes must be artificially introduced into the flow. Several methods of doing 
this successfully have been developed, and in all cases the change in density was 
made to coincide with the surface of discontinuity in velocity so that its effect on 
the flow may be assumed to be small. 

The first method employs a column of air heated by passing over a hot soldering 
iron. With the plate stationary the column was arranged to pass just outboard of 
the edge of the plate and was then barely discernible in the shadowgraph 
pictures. When the plate is moving, the heated air is drawn behind the lee side of 
the plate to be on the inner side of the surface of discontinuity, which is then 
shown up with great clarity. 
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Other methods employ, on the surface of the plate, fluids which evaporate. 
A plate of balsa wood was used and a narrow band of benzene applied by paint 
brush to the upper and lower surfaces. The benzene soaks into the balsa wood 
and then starts to vaporize. When the plate accelerates, air flowing over the 
surface increases vaporization of the benzene and the vapour is drawn into 
the vortex sheet. Because of the difference in density between the vapour 
and the air, the paths of their boundaries are clearly seen in the shadowgraph 
pictures. 

When the benzene is brushed on in a narrow band across the centre of the plate, 
normal to the light path, the resulting photographs have the effect of showing up 
the flow in a plane until the vortex and vapour sheet has moved far enough from 
the plate to be affected by end effects of the plate. 


3. Results 


The plates tested were all of 3in. square planform and each had a cross-section 
which was constant in the direction of the light path. A variety of plates was used, 
these differing in their cross-sections as may be seen in the photographs, 
figure 6(a)-(e). Also visible in each figure is a scale of pins at } in. intervals, the 
initial position of the plate tip being opposite the bottom pin. This, together with 
the time intervals indicated, defines the motion of the plate. Care was taken to 
ensure that the undulations of the vortex sheet were not caused by vibrations 
of the model by testing plates of different stiffness; nor were they caused by any 
three-dimensional periodicity across the span of the model. The absence of such 
periodicity was verified by simultaneously taking pictures in two planes at right 
angles to one another. 

The main sequence of pictures in figure 2 shows the growth of the edge vortex 
sheets behind a plate with a convex front surface and a flat rear surface as it 
accelerates to, and remains at, a speed of 24 ft./sec. Figure 3 shows the distance 
from rest of the plate as a function of the time taken. Figure 4 shows an enlarged 
section at ¢ = 5-41 msec. 

The sequence of photographs in figure 2 reveal the growth and the structure of 
the vortex sheets very clearly, and the small-scale formation of a large number of 
centres of vorticity within the big-scale rolled-up vortex sheet, in the manner 
described by Prandtl, is distinctly visible. The regularity of the process is note- 
worthy. Figure 3 also shows the movement of the core of the big-scale rolled-up 
vortex sheet. It is found that the velocity of the core is about half that of the 
plate. This is consistent with what is known about the corresponding three- 
dimensional flow past a rectangular plate of small aspect ratio where the upper 
edges of the tip vortex sheets are inclined to the mainstream at about half the 
angle of incidence (Bollay 1937, 1939; Mangler 1939; Kiichemann 1955). This was 
also found to be the case when the plate was accelerated to half the speed, i.e. 
12 ft./sec, although the initial undulations in the shear layer appear less pro- 
nounced as shown in figure 5. Future investigations may relate the change in the 
initial undulations to effects depending on the Reynolds number, which, at 
present, are not known. 

The remaining photographs, figure 6, show the pronounced effect of the plate 
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cross-sectional shape on the regularity of the distortion of the shear layer. + The 
differences in the flow are of secondary nature: a round edge, figure 6 (c), produces 
a less orderly flow with, possibly, two separation lines. A very pronounced 
shoulder, figure 6(e), leads to a secondary separation. Sharp edges which fix 
separation lines produce clean flows. 


| Plate velocity 24:0 fi./sec 


Inches 


» 
Vortex centre velocity 
——+10-9 ft./sec 
| | 
| | 


8 10 


Time (msec) 


FicurRE 3. Movement of the plate and vortex centre from rest 
obtained from figure 2(a)-()). 


4. Conclusions 

A technique has been presented which should prove useful in the examination 
of the development of shear layers and other time-dependent flows at low speeds. 
Results obtained so far clearly show the distortion of a shear layer into regular 
vortex filaments and finally into turbulent motion. Further experiments are 
necessary to understand the dependence on the Reynolds number of the regular 
manner in which the shear layer breaks up. 
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+t In some of the pictures, particularly figure 6(b), some particles of balsa dust are 
visible; these have been shed by the plate during the acceleration phase. 
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(a) ¢ = 1-05 msec; instantaneous velocity (v,) = 


(d) t = 4:30 msec; v, = 21-0 ft./see. 
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(7) t = 10-66 msec; v, = 24-0 ft./sec. 


FicurRE 2. Nequence of shadowgraphs of convex plate section with final velocity 
24 ft./sec. 
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Ficure 4. Enlargement of figure 2(e), ¢ = 5:41 msec, 
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Figure 5. Shadowgraphs of convex plate section with final velocity 12-1 ft./see. 
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Figure 6. Shadowgraphs of various plate cross-sections at ¢ = 6-53 msec. All 
velocities 24 ft./see approximately. 
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A turbulent equatorial jet 


By ROBERT R. LONG 
Department of Mechanics, The Johns Hopkins University, Baltimore 


(Received 31 January 1961 and in revised form 4 April 1961) 


A recent study of a laminar jet in a rotating spherical shell of fluid is extended 
to the case of a turbulent planetary jet at the equator of a rotating, stratified 
atmosphere or ocean. General forms of the velocity, density and pressure func- 
tions of both the mean motion and the turbulence are derived by a dimensional 
analysis applied to the mean and perturbation equations. The horizontal and 
vertical dimensions are estimated, based on the three characteristic constants 
of the problem, which are the momentum transfer, the stability and a rotation 
parameter. The estimates are in good agreement with the dimensions of the 
Cromwell current, i.e. the equatorial undercurrent of the Pacific Ocean. 

To the first order of approximation, the mean axial velocity in the theory is 
independent of distance along the jet axis. The mean horizontal transverse 
velocity component is much smaller and decreases upstream. The mean vertical 
velocity is extremely small, also decreasing upstream. The two horizontal velo- 
city components of the turbulence are of the same order and, in the undercurrent, 
are about one-fourth the mean axial velocity. The vertical turbulent component 
is much smaller. Finally, it is shown that the eddy-viscosity concept is inappro- 
priate for this problem because at least one of the eddy coefficients would have 
to be negative. 


1. Introduction 

In a recent paper (Long 1960) the author investigated a laminar, two-dimen- 
sional jet in a rotating spherical shell of homogeneous, viscous fluid. The flow is 
predominantly zonal, or west-east, and a suggested physical mechanism is a 
balance of vorticity brought into the jet by the weak north-south motion and 
diffused by friction. The order of magnitude of the horizontal width of the 
theoretical jet is given by 7/f62 ~ 1, where 4, is the width, w is a representative 
velocity, and / is the north-south rate of change of the Coriolis parameter. This 
expression was used to estimate the width of the Cromwell current or equatorial 
undercurrent. This is a west-to-east current at the equator in the Pacific Ocean, 
located a hundred metres or so below the surface (Cromwell, Montgomery & 
Stroup 1954; Fofonoff & Montgomery 1955; Hidaka & Nagata 1958; Knauss & 
King 1958; see also the entire volume Deep Sea Research, 6, 1960). If we use values 
u = 15x 10? emsec"!, = 2:6 x appropriate to the Cromwell 
current, we obtain 6, ~ 2:6 x 10’cm. This agrees well with the observed width, 
300 km, as reported by Knauss (1960). 

‘he laminar theory had two deficiences in relation to such geophysical jet 
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flows as the undercurrent and the equatorial atmospheric jet (Lettau 1956): 
(1) the theory assumed that all quantities were independent of height, whereas the 
natural jets have much smaller thicknesses in the vertical than in the horizontal; 
(2) the motion was laminar, whereas the natural jets are no doubt fully turbulent. 
It was argued that the results could be applied to the undercurrent by introducing 
the concept of eddy viscosity, but the validity of this concept is always doubtful. 

In this paper we investigate a fully turbulent jet in an infinite atmosphere or 
ocean. The axis of the jet is along the equator, and a general view is as pictured 
in figure 1. We may suppose that the flow is into a sink in the vicinity of x = 0, 
y= 0,z=0. We confine attention, however, to the flow in the vicinity of 
(— Xp, 0, 0) far upstream.t 


North Up 
y fe 
| g (24) 
u 
Equator — — Ss Equator 
South Down 


FicureE 1. Schematic picture of mean motion in jet. The transverse velocity fields may 
have signs different from those suggested by this drawing. 


2. Basic equations 
Weconsider the fluid to be incompressible but inhomogeneous. At great distances 
vertically and horizontally from the jet axis, the gradient of density is linear in 

the vertical; i.e. 


Oz g° 


> 0, 


where / is a constant, g is gravity, and p, is a reference density. We also make the 
Boussinesq approximation (Boussinesq 1903) that density variations are negli- 
gible in their effects on the inertia of a parcel and in the equation of continuity, 
but of fundamental importance when they multiply gravity. Then, if we define 
a perturbation density p* by 


kz 
P= Pot 


p* 2 
and G, P by G=-g—, P=—+92- ke, (1) 
Po Po 
the equations of the problem can be written 
vu, +wu,—fyv = —P,, (2) 
Uv, + vv, + wr, + = —P,, (3) 
uw, +ww,—-G = —P., (4) 
G,+uG,+vG, +uG,+kw = 0, (5) 
U,+v,+w, = 0. (6) 


t We consider only westerly jets. Qualitative arguments, attributed to C. C. Lin, 
indicate that westerly jets are stable, easterly jets unstable to large-wavelength per- 
turbations (Veronis 1960). 
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Equations (2) to (4) are the equations of motion, equation (5) expresses the in- 
compressibility assumption, and (6) is the equation of continuity. 

We have neglected molecular friction and diffusion. We know, of course, that 
these terms are important in so far as they affect the very small eddies, but it is 
well established that in fully developed turbulent flow at high Reynolds numbers, 
the effect on the mean motion as well as on the energy-containing eddies is 
negligible (Townsend 1956). 

The terms — fyv and fyu are the Coriolis accelerations, in which the Coriolis 
parameter 2Q sin @ (Q is the angular rotation speed of the sphere, # the latitude 
angle) has been approximated by 2Qy/a = /y, where a is the radius of the sphere. 
This is a very good approximation if the north-south extent of the jet is sufti- 
ciently limited. Indeed, the error is of the same order as that involved in using 
equations in Cartesian instead of spherical co-ordinates, namely o(A-*), where 
A is the ratio of the east-west and north-south length scales (Long 1960). 

Finally we notice that in the two-dimensional theory (Long 1960), the solution 
described a jet motion that could exist at any latitude. In this paper the jet 
must be located at the equator if the equations are to have the form of (2)-(6). 

We now take ensemble averages. Then, for example, we have 


(7) 
where % is the average velocity and uw’ = 0. If we substitute into (2)-(6), we get 
a set of averaged and unaveraged differential equations: 


uu, + vu, + wu, — = —(w'v’), — (w'w’)., (8) 
UD, + V0, + we, + = —P,—(u'v’), —(v'v’), — (9) 
UW, + vw, +ww,—G = (10) 
A-i An 
UG, + 04, + wG,+ wk = — (w’G@’)., (11) 
A=3 
U,+0,+ Ww, = 9, (12) 
and 
(WW), + +0U,+0U, — (Ur), 
A-t 
+w'u,+ wu, +w'u,—(u'w’).— (13) 
+ Uv, + — (w'v’), + + bY, —(v'v’), +0'v, 
Ash 
+ We, + — + (14) 


efA-? 
+w'w,+ vw, (15) 
30-2 


~ 
56): 
the 
tal; 
nt. | 
Ing 
ful. 
or 
red) 
of 
| 
| 
| 
) 
| 
n 
e 
. 
1 
e = 
) 
} 
| 
) | 
| 
| 
‘ 
Bi, 


468 Robert R. Long 
A-t 


+0'G, + wG,+ w'G, — (w'G’), + w'G,+ kw’ = 0, (16) 
U,+0,+w, = 0. (17) 


The symbols A-1, e?A-1, ete., over the terms in these equations express the order 
of magnitude of the terms in a way which will be explained later. 

If we now integrate equation (8) over a section of the jet, and assume that the 
velocities go to zero sufficiently fast at great distances from the jet axis, we get 
that 


Oo _ AS 
J= | | (P+wv+u'w' dydz (18) 


is constant, where we have used integration by parts and the equation of con- 
tinuity to transform the Coriolis term. The constant J is fundamental to the 
problem. It is called the momentum transfer. If we specified the conditions at 
the sink precisely and all other conditions of the problem, and found the exact 
solution corresponding to zero motion at |y|—> 00, |z| > 0, the integral (18) 
would be constant (independent of x), and J would be automatically determined. 
We confine attention, however, to the flow at large distances from the sink, and 
the solution we seek is not valid in the region (near the origin) where conditions 
are imposed that involve the precise description of the sink, say the flux Q and 
a linear dimension a. In adopting (18) as an alternative condition we must recog- 
nize the possibility that this may not be an adequate substitute for the two 
conditions involving a, Q, and that an indeterminacy may result. (Additional 
remarks about this point are contained below.) 


3. General form of solution 

We first make boundary-layer type estimates of the size of the terms in the 
mean equations. For example, it is plausible that the Coriolis and pressure terms 
alone dominate in (9), and that the vertical pressure variation is hydrostatic 
in (10). A first approximation for the unknowns in the equations of mean motion 
can then be obtained by using a generalized dimensional analysis, i.e. a search 
for invariance under a general affine transformation (Birkhoff 1950). We may 
then use the equations defining the Reynolds stresses to estimate the magnitude 
of such turbulent quantities as w’v’. A combination of these estimates of mean 
and turbulent quantities leads us to the general form of a first approximation 
to the solution and to a set of simplified differential equations that determine 
the first approximation. This in turn suggests a consistent scheme of successive 
approximations to the exact solution of the fully turbulent jet problem in the 
form of power series in certain small variables and parameters. We do not, to 
be sure, show that other expansions are not possible. 

We first shift the origin of our co-ordinate system from the orifice to an origin 
at a great distance toward the west, X,. Then 2 is the variable distance along 
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the x-axis from the section in question. We introduce a typical velocity and 


length = (19) 
L= piJikis, (20) 
and non-dimensional independent variables 
E=a/L, y=y/L, C=2/Le, 7 =telL, (21) 
where e = 22) 
We also need a non-dimensional function of x: 
A= (23) 


We assume that A is very large. It is of the order of the ratio of the axial scale of 
the mean motion to the width of the jet if we assume 7 ~ 1. We see that if the 
quantities € and 7 in (21) are both of the order of one, € is the ratio of the vertical 
and north-south length scales. We assume that ¢ is very small. (It is of the order 
of 10-3 in the case of the equatorial undercurrent.) 

We now search for solutions for the mean and turbulent quantities in the forms{ 


ule = U,+ +... +2U,4+..., (24) 

tAle +... (25) 
wAlec = +eW,+..., (26) 
Pie? = (27) 
= +..., (28) 
ww’ Aj? =F, (29) 
u'v' Ale? = N,+A-3N, +... +e2Ng+..., (30) 
u'w' = H,+ (31) 
Ale = (32) 
v'w'A/e’e = (33) 
w'w'A|c%e? = +..., (34) 
wG'eLA/c? = (35) 
= B,+ AB, +... +B (36) 


+ Although we define A with the idea in mind that X, represents the distance to the 
orifice from the section at which we perform our analysis, in fact nothing in the analysis 
will be changed if we regard X, simply as a sufficiently large positive quantity with the 
dimensions of length. Some consequences of this are discussed below. 

t The first approximations to the solution represent self-preserving flow, but in the 
development of this paper we show that this type of flow is possible at great distances from 
the orifice by considering all equations and conditions of the problem, not simply the 
equations of mean motion. Since the analysis may be applied to jets and wakes in homo- 
geneous, non-rotating fluid, the approach of this paper removes the necessity for the 
customary assumption of self-preserving flow. The point is discussed by Townsend (1956). 
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= C,+A-4C,4+.. 
w At/e = u,+A-tu,+ 
Atle = v, + + 


w'At/ce = w, + +. 


= 


G'eLAt/c2 = g, +. 


+. 


+...; 


Wey +..., 


+..., 


+675) +..., 


where the functions denoted by capital letters, U,, 3, ete., are functions only 
of 7, ¢, and those denoted by small letters, w,, v,, etc., are functions of £, 9, ¢, 7. 

If these expressions are substituted into equations (8)—(18), we find that the 
sets of quantities U,, Vi-p,, g,; U>, Vi-ps, Jo; etc., satisfy sets of equations in 
independent variables , 7, ¢, 7, from which all the constants of the problem are 
missing and which do not involve A or e. If we assume, therefore, that U,, V—p,, 9, 
as well as €, 9, ¢, 7 are all of order one, we may then estimate the magnitude of 
all terms in equations (8)—(18). A quantity such as A-}, e?A-*, ete., written over 
a term is the ratio of that term to the dominant terms of the equation. 

The first set of equations is 


= (44) 

1={" (?, dt, (48) 
= = (50) 
2, (51) 

Bee +w, = 0, (52) 
= 0. (53) 
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These are 10 differential equations in 14 unknowns U,, Pi, Gy, 
Pris 91» Ny, Hy, By, Cy. We can make the system determinate, however, by adding 
the four equations 


a 

N, = lim u,v, dé, (54) 
A> «0 A 0 

12 

H, = lim uw, dé, (55) 
A> 2 A 0 
a 

B, = lim = | 1491 4£, (56) 
A> A 0 
1 A 

C, = lim | dé. (57) 
A> A 0 


Since .V,, H,, B,, C, do not vary with £, we substitute for the ensemble average 
a space average along the z-axis. The equations (43)—(47), (49)—(53) and (54)—(57) 
are now a complete set of equations. Notice that if we multiply (52) by g, and 
average, we get 

B, i) =0. (58) 
This may be substituted for either (56) or (57). 

The mathematical problem of the first approximation is not completely stated 
by the system of equations (43)-(57). Indeed, we see that if U,, V,, W, Py, G4, 
Uy, U1, Wy, Py; Jr, Ny, Hy, By, C, is a solution of the system, then there is an infinity 
of solutions U,, «7V,, a?W,, P,, Gy, av, aw, ag,, «7H, B,, 
where « is any constant. Notice that as « is varied the solutions change in just 
the same way as a given solution changes when X} is varied. (The difference 
between X,—2 and X, in A is negligible for the purposes of the first approxi- 
mation.) A similar indeterminacy arises in the non-rotating, homogeneous, 
turbulent jet. If the eddy-viscosity concept is used, it appears as an indeter- 
minacy in the eddy Reynolds number (Townsend 1956). 

One reason for the indeterminacy is that we have said nothing yet about our 
origin of time or about the turbulent velocity, pressure and density fields that 
exist at that instant. Something must be specified about those turbulent fields, 
but since we have assumed no variation with time of the mean quantities the 
initial fields are not completely arbitrary. For example, if we assume solutions 
for the turbulent quantities in the form 


Uy = Uy, cos (a& + br) + uy, sin (af +57), 


01, sin (af + b7) + v4, cos (af + b7), 


v. 1 


= W,, sin (af +b7) + wy, cos (ag + br), 


| 


Py = Pi, 008 (a& + br) + sin (ag + 


91 = 911 008 + br) + 942 8in (af + b7), 
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where U1, Uj, etc., are functions of 9, €, and a, b are constants, we find from 
(54) to (57) . 

in 

“wi 
Ch = 
On substitution into (43)-(53), the sine and cosine terms cancel out and we get 
15 differential equations in 15 unknowns Uj, Vj, ..., In 
this way we can arrive, in principle, at a solution of the first-approximation 
problem, but theinitial conditions that are implied by this solution are very special. 
It is possible, of course, to assume solutions for the turbulent quantities that will 
satisfy more general initial conditions, but this increases the complexity of the 
equations that must be solved. 


4. General properties of the jet 

If we suppose that all the non-dimensional quantities that we have introduced 
(except A, €) are of the order of 1, and if we confine attention to the first approxi- 
mation in (24)—(42), we can state a number of properties of the jet: 

(1) Combining U, ~ land 7 ~ 1, we find that 6, ~ «3-4, where d, is the width 
of the jet. The width is small if the speed is low and if the rotation is high. In 
order of magnitude this is the same as the width of the laminar jet (Long 1960), 
although the laminar jet widens gradually with distance from the sink, whereas 
the width of the turbulent jet is uniform. 

(2) Combining £ ~ 1, 7 ~ 1, we find that 6, ~ wk-*, where 6. is the vertical 
thickness of the jet. Hence the vertical extent of the jet is small if the speed is 
low and the stability large. The vertical thickness is also independent of distance 
from the sink. Notice that the motion so adjusts itself that the non-dimensional 
numbers representing the rotation effect and the stability are of the order of one. 

(3) The axial mean velocity is constant along the length of the jet. 

(4) The transverse mean velocities both decrease upstream as 1/X,. In magni- 
tude, < wand w < v. 

(5) The mean-density field is independent of distance along the axis. 

(6) The two horizontal turbulent velocity components are of the same order 
and smaller than the axial mean velocity. The vertical and horizontal transverse 
turbulent velocities are larger than the corresponding mean quantities. The 
vertical turbulent velocity is much smaller than the horizontal components. 

(7) All turbulent quantities decrease rather gradually upstream as Xj?. 

(8) The horizontal scale of the turbulence is of the order of the horizontal 
width of the jet; the vertical scale of the turbulence is of the order of the thickness 
of the jet. 

An important result may be obtained concerning the eddy-viscosity concept 
applied to this problem. For example, if we assume the turbulent transports are 
proportional to the mean density gradients 
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where K,,, AK, are horizontal and vertical components of eddy diffusivity, the 
form of (46) is preserved only if 


1 1cG 
1 1 (eG 


where S,,, S, are non-dimensional functions of 7, ¢. The coefficients of diffusion 
therefore decrease upstream as we would expect, and in general vary with the 
constants of the problem in a plausible way. But using (58) we have 


1 1 (0G, .\2 


This can only be true, however, if the diffusion coefficients are allowed to be 
negative. We see clearly that an analogy with molecular diffusion and friction 
is not useful in this problem of turbulence. 


5. Circumferential equatorial jets 

In the next section we apply the theory of the turbulent jet to the equatorial 
undercurrent. Before doing this we may mention here that westerly currents 
at the equator seem to be a very common feature of the atmospheres of the 
planets and the sun. In astrophysics it is known as the equatorial acceleration 
(Wasiutynski 1946; Cowling 1953; Hess 1951). In Jupiter, for example, this 
acceleration occurs in a band between —7° and +7° latitude and amounts to a 
relative angular velocity of about 1°, of the angular velocity of the planet 
(w/Q ~ 0-01). In the atmosphere of the sun the acceleration is very strong, but 
very variable within a sunspot cycle, @/2 varying from 0-3 to 0-7 or so. The band 
containing high-speed fluid is between — 55° and + 55° when w/Q is a minimum; 
when @/Q is a maximum, the width of the band is very indefinite, but certainly 
even broader than this. 

The same phenomenon occurs in our own atmosphere, although it is not yet 
known whether it extends round the entire globe (Lettau 1956). This jet is called 
the Berson westerlies by Lettau in honour of Berson who first reported the 
phenomenon (1910). Observations indicate that it is located in the stratosphere 
at around 20km. It isabout 5km thick and 1500 km in north-south extent. The 
speeds are of the order of 10% cm sec™. 

The theory of this paper can be modified to apply to this phenomenon. We 
need only recognize that there can be no variation of the mean quantities with 
longitude, so that A is simply an undetermined large number. The order of the 
width of the jet is the same, 6,, ~ utp-4, or 60 ~ (w/2Q)!. This seems to describe 
in a rough way the variation in the size of the jet in the astrophysical examples. 

If we apply the theory to the Berson westerlies, we get 0, ~ 6 x 10° em. This is 
considerably less than the observed width, but our knowledge of this current is 
still quite limited. If we use wu = 108, k = 4x 10-4 we find 6, ~ 5 x 104em which 
is very much too small. However, we should note that the westerly jet is im- 
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bedded in a belt of easterlies with a general drift of perhaps 3 x 10? cm sec-!, 
It is more reasonable then to use for u the relative velocity of 4 x 10°. This yields 
6, ~ 2x 10°cm and 6,, ~ 1-2 x 108 em. 


6. Application to the equatorial undercurrent 


In any attempt to apply the above theory to the equatorial undercurrent we 
must recognize that an ocean surface exists in one case and not in the other. 
Our theory of a jet in an infinite fluid can be adapted to this case, however, as 
follows: if we guess that the mean theoretical axial velocity profile in a vertical 
section resembles that in figure 2 (counter currents are typical phenomena in 
rotating and stratified fluids), we could put in a free surface S at the indicated 
position without changing the situation substantially, provided the surface 
stress in the natural phenomenon is substituted for the stress exerted by the 
fluid above S in the theoretical situation. If we adopt this viewpoint, this force, 
supplied by the wind stress in the case of the ocean, may be related to the depth 
of the jet axis below the surface. 


FicureE 2. Schematic picture of vertical velocity profile of undercurrent. 


In applying our theory to the Cromwell current we take as observed 
values (Knauss 1960) a mean velocity of 1-5 x 10? em sec~!, width 3 x 107 cm, 
vertical thickness 2 x 104 cm, longitudinal length 5 x 10°cm or more, £ = 2:3 x 
10°3em-'sec!. If we use the given values of uw and /, we may compute 
6, and We get 6, ~ 2:6 x 10’cm, 6, ~ 104cm, where we have estimated the 
stability parametert k = (g/p)dp/dz ~ 4x 10-4sec~*. Since the estimates are 
only supposed to be order of magnitude estimates, these results are remarkably 
close to the mark. The value of A is about 20 in this case, so that all turbulent 
quantities are } to } of the corresponding mean quantities. 

Although we do not solve for the mean velocity and temperature fields in the 
jet we may make certain inferences about them. In figure 3 we picture the flow 


+ This estimate was supplied to me by Professor R. B. Montgomery. 
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in a vertical north-south cross-section. Since the motion is westerly in the jet 
the y- and z-equations of mean motion (or the thermal wind equation) require 
a temperature field in which the isotherms are roughly as drawn. This suggests 
the vertical and horizontal transverse velocities shown, and this is consistent 
with the equation of continuity (47). It is important to remark that this is 
precisely the observed form of the constant density surfaces in the vicinity of 
the Cromwell current as we see in figure 4. This transverse velocity field is also 
consistent with the Ekman wind-drift theory, which requires that the surface 
water move north and south away from the equator in this region. 
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Figure 3. Mean temperature distribution in jet. 


The author wishes to thank Professor R. B. Montgomery and Mr R. 8. Arthur 
for a number of stimulating conversations about the undercurrent. This research 
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Appendix 

Jets are very common in stratified and rotating fluid systems. For example if 
a rotating fluid is heated, very pronounced, irregular jets occur at the free surface 
(Fultz et al. 1959). They resemble the famous Jet Stream of the atmosphere and 
other currents like the Gulf Stream in the oceans. In these phenomena the 
rotation itself, rather than the variation of the rotation, is probably the more 
important effect. An analysis, essentially the same as that above, can be applied 
to such a jet, using the same equations as (2)-(6) except that f = 2Q replaces /y. 
The solutions all have the same form as (24)-(42). However, now we have 
c= Jifiks, L = Jif-ikt, ¢ = fk-?, and in the differential equations of the first 
approximation 7 is missing in the terms corresponding to the Coriolis forces. 
We may compare these results with observations of the Jet Stream and Gulf 
Stream (Newton 1959; Stommel 1960). For example, Newton’s data give for 
the ratio of depth to width 0-005 for both. We get precisely 0-005 for e€ if we use 
the reasonable values f = 10-4, k = 4x 10-4. On the other hand, the theoretical 
estimates of the dimensions of the Jet Stream and Gulf Stream are too low by 
a factor of 2 or so. 
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FiGuRE 4. Specific-volume anomalies in centilitres per ton on vertical section in Pacific 
Ocean at 150° W. in July-August 1952. Observations from Hugh M. Smith of Pacific 
Ocean Fishery Investigations, U.S. Fish and Wildlife Service. Drawing from unpublished 
manuscript by R. B. Montgomery and E. D. Stroup. 


REFERENCES 


Berson, A. 1910 Bericht iiber die aerologische Expedition des Kénigl. Aeronautischen 
Observatorium nach Ostafrika im Jahre 1910. Ergebn. Arb. Preuz. Aeronaut. Obs. 
Lindenberg: Braunschweig. 

Brrxuorr, G. 1950 Hydrodynamics, chapter 1v. Princeton University Press. 

BovusstnesqQ, J. 1903 Théorie analytique de la chaleur, Vol. 2, p. 172. Gauthier- Villars. 

CowLinG, T. G. 1953 Solar electrodynamics. Article in The Sun (Ed. G. P. Kuiper). 
University of Chicago Press. 

CROMWELL, T., MonTGOMERY, R. B. & Stroup, E. D. 1954 Equatorial undercurrent in 
Pacific Ocean revealed by new methods. Science, 119, 648. 


F 
| e | \ | \ 
| 
e e \ 
e e e 
e 
W 


en 


A turbulent equatorial jet 477 


Foronorr, N. P. & Montcomery, R. B. 1955 The equatorial undercurrent in the light 
of the vorticity equation. Tellus, 4, 518. 

Futtz, D., Lone, R. R., OwEns, G. V., BoHan, W., Kaytor, R. & WEIL, J. 1959 Studies 
of thermal convection in a rotating cylinder with some implications for large-scale 
atmospheric motions. Meteor. Monograph, 4. 

Hess, 8. L. 1951 The atmospheres of the other planets. The Compendium of Meteorology, 
p- 391. 4 

Hipaka, K. & Nacata, Y. 1958 Dynamical computation of the equatorial current system 
of the Pacific, with special application to the equatorial undercurrent. Geophys. J. 
Roy. Astr. Soc. 1, 198. 

Knauss, J. A. 1960 Measurements of the Cromwell Current. Deep Sea Res. 6, 265. 

Kwauss, J. A. & Kine, J. F. 1958 Observations of the Pacific equatorial undercurrent. 
Nature, 182, 601. 

Lettau, H. 1956 Theoretical notes on the dynamics of the equatorial atmosphere. 
Beitr. Phys. Atm. 29, 107. 

Lone, R. R. 1960 A laminar planetary jet. J. Fluid Mech. 7, 632. 

Newton, C. W. 1959 Synoptic comparison of Jet Stream and Gulf Stream systems. 
The Atmosphere and Sea in Motion, p. 288. New York: The Rockefeller Institute 
Press. 

StoMMEL, H. 1960 The Gulf Stream, p. 202. University of California Press. 

TownsEnD, A. A. 1956 The Structure of Turbulent Shear Flow, p. 89. Cambridge Uni- 
versity Press. 


VeRONIS, G. 1960 An approximate analysis of the equatorial undercurrent. Deep Sea 


Res. 6, 318. 
WasrutyNskr, J. 1946 Hydrodynamics and structure of stars and planets. A strophysica 
Norvegica, 4, 497. 


| 
| 
| 
| 
| 
| 
| 
fic 
fic 
ed 
IS. 
in 


478 


REVIEWS 


International Summer Course in Plasma Physics. Danish Atomic Energy 
Commission, 1960. 645 pp. Dan.kr. 60,00. 


In August 1960, the Danish Atomic Energy Commission organized a two-week 
course on plasma physics at their research establishment Ris6. A distinguished 
team of a dozen lecturers was assembled from half a dozen countries, and over 
130 students were enrolled. Only three months later the text of the lectures was 
published, in a report containing about 200,000 words. Evidently there is much 
interest in plasma physics these days, and there also seems to be a great deal to 
say about it. A brief reference, at the end of the report, lists the titles of fifteen 
additional seminars which took place, but whose text is not reproduced. Surely 
no group of students or lecturers can ever have worked harder. 

Now the study of plasmas is interesting, for many reasons. At first sight a 
plasma may appear to resemble a gas, for it is compressible and has no structure. 
But the forces between the charged particles in a plasma are electromagnetic, 
and act over longer ranges than the forces between the molecules in a gas. 
Therefore we can normally consider that molecules interact only intermittently, 
when they approach close enough to each other, and that then they interact 
strongly. But it is more profitable to think of the particles in a plasma as inter- 
acting with each other the whole time, though in most cases the interaction 
between any pair of plasma particles is weak. Here lies the basic difference 
between these two states of matter: where close collisions dominate, the out- 
come of a particular encounter is rather uncertain, and this introduces a random- 
ness into the mechanies of a gas, which is less pronounced in a plasma. A plasma 
particle interacts, essentially, with a continuum, whose behaviour is more 
predictable than that of a single molecule. We can therefore follow its progress 
in phase space more easily. A plasma thus tends to be more internally coherent, 
and a whole spectrum of waves becomes possible which is absent in a gas. The 
waves are important because they can lead to the growth of non-uniform electric 
and magnetic fields. 

A plasma will respond strongly to an applied electric or magnetic field, and 
will at first try to exclude it from its interior, though in the case of a magnetic 
field this often leads to unstable configurations. In a plasma laced with a 
magnetic field the paths of the charged particles twist around the lines of force. 
If high-energy electrons are present, and they seem readily to be produced in 
such an active medium, they become powerful emitters of electromagnetic 
radiation. 

Even the process of turning a gas into a plasma raises intriguing problems. 
Though they may have energy enough to detach electrons from one another, 
atoms and molecules are inefficient in producing ionization. But once a few free 
electrons exist and have picked up enough energy, the ionization sweeps on 
like an avalanche. Other interesting effects are the process by which neutrals 
and ions can exchange electrons, and the way in which contaminants, particularly 
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ions with one or three electrons, can be excited and then radiate. These are both 
important mechanisms for cooling a plasma. 

Plasmas can be studied in many and widely different physical situations. They 
occur, for example, in radio stars and in the tails of comets, in the Van Allen 
belts and the curtains of the polar aurora, in lightning flashes and in television 
tubes. Yet none of these applications provoked the intensive modern work on 
plasmas. This is due to the hope that somehow one might be able to confine for 
a long enough time a plasma which is sufficiently hot, say at a temperature of 
about 108 degrees K. One could then set up a controlled thermonuclear reaction 
(c.T.R.), make helium from deuterium or tritium, and thus release a vast amount 
of energy. The confinement must be achieved by a skilfully designed magnetic 
field. The background to the Ris6 lectures is the question whether this can be 
done. 

Here must be the hardest technical problem that man has yet set for himself. 
It is worth thinking about, for after its solution our fuel requirements will be met 
for the indefinite future. Though no one can reasonably maintain that it is an 
urgent question, it was, in the curious political circumstances of the last decade, 
turned into the theme of an international competition, rather like the race 
towards the moon. Fortunately it has now been realized that the fundamental 
physics must first be understood. This is quite well reflected in the various 
lectures. The exception is one contribution, describing a thermonuclear reactor, 
whose author first disarms his potential critics by admitting that plasmas are 
too unstable for his liking. His machine is therefore so designed that the plasma 
is not given time to perform its tricks properly, and by hurrying to forestall it, 
his machine might just scrape by with a very slight profit. It is all rather touch- 
and-go, like the old problem of filling a leaky bath tub with both taps full on. 

In general the lecturers did not take such a point of view. They may have been 
thinking about thermonuclear reactors, but in doing so they felt obliged to— 
consider almost all the different aspects of the behaviour of a plasma. 

The contributors seem to have had different ideas of their responsibility to 
the audience. Some are represented in the report by photographic copies of 
papers published in well-known and highly accessible journals, none of them less 
than two years old at the time. Should these have been reproduced in the 
report? Surely anyone who wishes to read them can go to a library and do so. 
Surely also it is the duty of a lecturer to do more than read to his students 
extracts from the classics. He should produce the occasional enlightening 
example, point out the critical stages in an argument, be less insistent on rigour 
and completeness, and in general make the audience think. 

Others again reported on very special studies, recently completed in their 
laboratories. These may indeed be most significant, but the reader does not 
know why and the audience was probably not certain either. 

Such criticisms apply only to the minority, however. The majority of the 
lecturers seem to have given accounts of their subjects as proper and as balanced 
as our present knowledge allows. It is an engaging feature of this book that in 
some cases the reports reproduce the lectures almost as they must have been 
given, complete with the occasional slip in algebra or in spelling, and including 
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a few samples of typical lecture jokes. After reading them for a while one feels 
as though one were present oneself in the classroom at Riso. 

Some of the contributions cause one to feel the same excitement that the 
audience must have felt on hearing physical problems discussed so lucidly and 
so graphically. Thus in the first lecture Rosenbluth describes, very briefly but 
with much insight, the basic dynamical processes of a plasma, and evaluates the 
characteristic quantities one has to watch. This must be the best introduction 
anywhere to these topics, and should be compulsory reading for all future plasma 
physicists. 

This exceptionally high level of presentation is not reached again—after all, 
how could it be? But Rosenbluth scores again in an excellent account of micro- 
instabilities; other first-class contributions include Thompson’s survey of 
instabilities, Simon’s discussion of the collisionless Boltzmann Equation, and 
Wharton’s three lectures on plasma waves and diagnostics. 

After digesting the various reports one is deeply grateful to the lecturers for 
all the new points of view and amazed at the number of physical problems they 
have raised. Then one begins to wonder whether the material is not too much 
biased towards the interests of the c.T.R. physicist. To think so would be an 
unfair criticism. Suppose that a dozen astrophysicists had organized the course. 
In essence they would have had to discuss the same types of problem, for the 
nature of the Coulomb interaction, the Boltzmann Equation, the processes of 
ionization and charge exchange are quite universal. The main difference is that 
c.T.R. work has sometimes tended to concentrate on the study of some very 
special configurations, and has thereby distorted the development of the subject. 
Let a c.T.R. physicist have the last word (Bickerton, on page 453). ‘The past 
experiments have been carried out in systems which were known or are now 
known to be unstable so that plasma behaviour has been complex and in detail 
incomprehensibie. In the future it is to be hoped that by shifting the complexity 
to the external apparatus and by making a serious attempt to set up stable 
conditions in the plasma the resulting phenomena will be simpler and more 


pleasing.’ F. D. Kann 
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THIRD INTERNATIONAL SYMPOSIUM ON 
RAREFIED GAS DYNAMICS 


This symposium will be held at the University of Paris, June 26-29, 1962, and 
will be similar in scope and purpose to the first (Nice, 1958) and second (Berke- 
ley, 1960) symposia with corresponding titles. 

The programme will range from topics of immediate significance for upper 
atmosphere and space flight to basic scientific studies and will inelude the 
following areas: 


Studies of the limits of the continuum theory or the quasi-equilibrium 
kinetic theory of gases 


Problems in kinetic theory of gases, particularly attempts to solve the 
Boltzmann equation 


Free-molecule and near-free-molecule flow in neutral and ionized gases 
The physies of surface interactions between gases and solids 


Boundary conditions for rarefied gas equations, and slip flow 


Experimental techniques and instrumentation developments bearing 
on the above, whether applied to laboratory or field experiments 


Edmond A. Brun (University of Paris) and Immanuel Estermann (Office of 
Naval Research, London) will act as co-chairmen: the editorial committee will 
be headed jointly by P. Germain, Paris, and L. Talbot, Berkeley. A call for 
papers will be distributed later. 


Inquiries from the U.S. should be addressed to L. Talbot, Department of 
Aeronautical Sciences, University of California, Berkeley, and from Europe to 
Laboratoire d’Aerothermique, 4c, route des Gardes, Meudon (S & O), France. 
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